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ABSTRACT 

The 7 -ray sky can be decomposed into individually detected sources, diffuse emis- 
sion attributed to the interactions of Galactic cosmic rays with gas and radiation fields, 
and a residual all-sky emission component commonly called the isotropic diffuse 7 -ray 
background (IGRB). The IGRB comprises all extragalactic emissions too faint or too 
diffuse to be resolved in a given survey, as well as any residual Galactic foregrounds that 
are approximately isotropic. The first IGRB measurement with the Large Area Tele- 
scope (LAT) on board the Fermi Gamma-ray Space Telescope (Fermi) used 10 months 
of sky-survey data and considered an energy range between 200 MeV and 100 GeV. Im- 
provements in event selection and characterization of cosmic-ray backgrounds, better 
understanding of the diffuse Galactic emission, and a longer data accumulation of 50 
months, allow for a refinement and extension of the IGRB measurement with the LAT, 
now covering the energy range from 100 MeV to 820 GeV. The IGRB spectrum shows a 
significant high-energy cutoff feature, and can be well described over nearly four decades 
in energy by a power law with exponential cutoff having a spectral index of 2.32 ± 0.02 
and a break energy of (279 ±52) GeV using our baseline diffuse Galactic emission model. 

The total intensity attributed to the IGRB is (7.2 ± 0.6) x 10 -6 cm -2 s -1 sr _1 above 
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100 MeV, with an additional +15%/— 30% systematic uncertainty due to the Galactic 
diffuse foregrounds. 

Subject headings: gamma rays: diffuse background, diffuse radiation 


1. Introduction 


The universe is filled with electromagnetic radiation, which can be characterized by a cosmo- 
logical energy density and spectrum. This extragalactic background light (EBL) is energetically 
dominated by thermal relic radiation from the last scattering surface observed as the cosmic mi- 
crowave background. Different physical processes characterize the EBL in each waveband — starlight 
in the optical, thermal dust emission in the infrared, and emission from active galactic nuclei (AGN) 
in X rays. The extragalactic 7 -ray background (EGB) provides a non-thermal perspective on the 
cosmos, which is also explored through the cosmic radio background as well as extragalactic cosmic 
rays (CRs) and neutrinos. 


The EGB represents a superposition of all 7 -ray sources, both individual and diffuse, from the 
edge of the Milky Way to the edge of the observable Universe, and is thus expected to encode di- 
verse phenomena (see Dernier (2007) for a comprehensive review). Guaranteed contributions arise 
from established extragalactic 7 -ray source classes including AGN, star-forming galaxies, and 7 -ray 
bursts. The beamed emission from blazars is sufficiently bright that statistically large samples of 
individual sources have now been detected to cosmological distances ( Ackermann et al.| | 2011 b). Ac- 
cordingly, the flux distribution of blazars even below the detection threshold for individual sources 
can in principle be estimated from a relatively firm empirical basis through an extrapolation of the 


observed flux distribution (e.g., 

Abdo et al. 

2010 c 

Ajello et al. 2012 , 2014), although a consensus 

has not yet been reached (e.g., 

Singal et al 

2012 

). For other populations, such as star- forming 

galaxies 

Pavlidou & Fields| 2002 j Thompson et al. 2007; Fields et al. 2010; Makiya et al. 2011| 

Stecker & Venters 2011 Ackermann et al. 2012c) a 

nd AGN with jets oriented obliquely to our line 
imulative intensity is almost entirely unresolved 

of sight ( 

Inoue 2011 Di Mauro et al. 2014), 

the cr 


by current instruments; calculations of the flux distribution incorporating physical models and/or 
multiwavelength scaling relations must be invoked to estimate their EGB contributions. There are 
additional theoretically well-motivated extragalactic source classes too faint to have been individu- 
ally detected thus far, including galaxy clusters and their associated large scale structure formation 
shocks ( Colafrancesco &; Blasi||1998 Loeb Sz Waxman]|2000| ) . 

At energies > 100 GeV, the interaction length for 7 rays with photons of the UV/optical/IR 
EBL becomes much shorter than a Hubble length, thus defining an effective 7 -ray horizon. The 
electromagnetic cascades initiated by both very-high energy 7 -rays (Coppi & Aharonian 1997) and 
ultra-high energy CRs (Berezinskii &; Smirnov 1975) interacting with the EBL create truly diffuse 
EGB contributions. 
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Finally, more exotic processes such as dark matter annihilation/decay may be present, though 
as yet unrecognized, in the EGB (Bergstrom et al. 2001; Ullio et al. 2002 [Taylor <fe Silk||2003 ). 

From an observational standpoint, there are two main challenges in measuring the EGB. One 
is to model the diffuse Galactic emission (DGE) created by CR interactions with interstellar gas 
(ISG) and interstellar radiation fields (ISRF), which is comparable to the EGB intensity at energies 
> 1 GeV even at the Galactic poles, and therefore represents a strong foreground to the EGB 
measurement. The second challenge is separating cosmic 7 rays from CR induced backgrounds at 
the detector level. For instruments in low Earth orbit, the CR intensity can exceed that of the 
EGB signal by a factor of up to ~ 10 6 . In addition there is a sizable flux of secondary particles 
that are produced by interactions of CR in Earth’s atmosphere. 


(Bergstrom et al. 2001 ; Ullio et al. 2002 [Taylor &; Silk||2003 


The existence of all-sky 7 -ray emission was first realized experimentally using 621 candidate 
7 rays collected by the OSO-3 satellite (Clark et al.|1968 Kraushaar et al. 1972), while Fichtel et al 


(1975, 1978) reported the first spectral measurement of an isotropic diffuse background using the 
SAS-2 satellite. Analyses using more sensitive instruments capable of detecting individual extra- 
galactic sources began reporting the residual all-sky average intensity after subtracting individual 
sources and DGE templates (e.g., Sreekumar et al. 1998; Strong et al. 2004, using EGRET). The 
remaining emission component is found to be approximately isotropic on large angular scales and 
is commonly called the isotropic diffuse 7 -ray background (IGRB). The sum of the IGRB and indi- 
vidually resolved extragalactic sources represents an upper limit to the total EGB intensity, since 
residual unresolved Galactic emissions may be present in the IGRB. For example, CR interactions 
with gas (Feldmann et al. 2013) or radiation fields (Keshet et al. 2004) in the extended halo of 


the Milky Way, unresolved Galactic sources such as millisecond pulsars (Faucher-Giguere & Loeb 


2010), as well as CR interactions with solar system debris (Moskalenko Sz Porter 2009) and the 


solar radiation field (Moskalenko et al. 2006; Orlando &; Strong 2007, 2008) have been considered 
as sources of approximately isotropic emission on large angular scales. 


The intensity attributed to the IGRB is observation- dependent because more sensitive instru- 
ments with deeper exposures can extract fainter extragalactic sources, whereas the total EGB 
intensity (assuming complete subtraction of all Galactic emissions) is the fundamental quantity. 


The Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (Fermi) is 
the first instrument with sufficient collecting area and CR-background rejection power to measure 
the IGRB at energies > 100 GeV. Since launch into low-Earth orbit on 11 June 2008, the LAT 
has operated primarily in a sky-survey mode that combined with a large field of view (2.4 sr) and 
good spatial resolution (~ 1° at 1 GeV) has enabled the most detailed studies of the DGE to date. 
The LAT is a pair-conversion telescope consisting of a precision tracker and imaging calorimeter, 
which are used together to reconstruct 7 -ray directions and energies, and a surrounding segmented 


anti-coincidence detector (ACD) to identify charged particles entering the instrument. Atwood 


et al. (2009) provide a description of the Fermi mission and LAT detector as well as details of the 


on-orbit calibration. Ground data processing, event selection, and instrument response functions 


are provided in Abdo et al. (2009c), Ackermann et al. (2012d), and Ackermann et al. (2012f). 

A first measurement of the IGRB spectrum between 200 MeV and 100 GeV based on 10 months 
of LAT data was published in Abdo et al. (2010b). In this paper, we present a refinement and 
extension of that analysis based on 50 months of sky-survey observations. Multiple improvements 
in event classification, Galactic foreground and CR background models, and analysis techniques 
have been implemented. Together with increased statistics, these updates allow for an extension 
of the LAT IGRB measurement by over a decade in energy, now covering the range from 100 MeV 
to 820 GeV. 


2. Data samples 


50 months of LAT data recorded between 5 August 2008 and 6 October 2012 are used for 
this analysis, corresponding to a total observation time of 1239 day:0 The events have been 
reprocessed with an updated instrument calibration which improves the agreement between data 
and simulation of the energy reconstruction quality, the point spread function (PSF), and certain 


classification variables, and thereby reduces systematic uncertainties (Bregeon et al. 


2013 


The LAT IGRB analysis poses especially stringent requirements on the 7 -ray purity of the 
event selection since both the signal and CR-background spatial distributions are quasi-isotropic. 
The residual CR background contamination must be reduced to a relatively small fraction of the 
total isotropic intensity in order to measure the IGRB with acceptable systematic uncertainty 
because the (not perfectly known) CR background is directly subtracted from the total isotropic 
intensity in the final step of evaluating the IGRB. 


The predefined event classes publicly available from the Fermi Science Support Center includ- 
ing P7ULTRACLEAN have insufficient CR background rejection performance for the IGRB analysis 
energies below E < 400 MeV and energies above E > 100 GeV. Therefore, we developed two dedi- 
cated event samples for the IGRB analysis with distinct selection criteria at low and high energies. 
The IGRB intensity measurements reported in Section [5] use the ‘low-energy’ sample for the energy 
range 100 MeV - 13 GeV, and the ‘high-energy’ sample for the energy range 13 GeV - 820 GeA^} 
Multiple event classifications are necessary in order to obtain the best possible compromise be- 
tween statistics and low CR backgrounds across the full LAT energy range since the compositions 
and interactions of CR backgrounds in the low- and high-energy regimes are rather different. The 
modifications to the baseline P7ULTRACLEAN classification for the two event samples used in this 
work are described below, and summarized in Table [l] 


RAT data recording is disabled for « 13% of the total on-orbit time, during passages through the South Atlantic 
Anomaly (SAA), a region with extremely high charged particle backgrounds. Only observation periods that passed 
data quality monitoring and where the angle between the LAT z-axis and the Zenith was below 52° are used for 
this analysis. Note that the actual livetime is ~ 8% smaller than the 1239 days observation time quoted here due to 
instrumental deadtime associated with event latching and readout. 

2 The reprocessed data are available from the Fermi Science Support Center (http://fermi.gsfc.nasa.gov/ssc/) 
together with a list of caveats regarding their usage. 

3 Each sample includes events from the full LAT energy range. The labels ‘low-energy’ and ‘high-energy’ refer 
to the energy regime for which the event classifications have been optimized. The energy overlap between samples 
allows for consistency checking between the low-energy and high-energy analyses (described in Section [ 3 J, a feature 
we used to verify that the specific choice of crossover energy between 5 GeV and 50 GeV does not affect the accuracy 
of our quoted results for the IGRB intensity. 
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2.1. Event selection 


The low-energy sample is a strict subset of events classified as photons according to the 


P7ULTRACLEAN event class definition (Ackermann et al. 2012f). To reduce the residual background 


of secondary electrons, positrons, and protons produced by CR interactions in the Earth’s atmo- 
sphere, which are the primary concern in the low-energy IGRB analysis, the following additional 
criteria are imposed: 


Tracker veto I: Part of the tracker is used as an additional veto to complement the ACD. Specif- 
ically, we require the reconstructed 7 -ray trajectory to cross at least two layers of active 
silicon strip detector area without producing hits in these detectors. This selection criterion 
significantly increases the efficiency of vetoing charged particles entering the LAT. 

Tracker veto II: We discard events for which the reconstructed pair-conversion vertex lies in 
the three x-y double-layers of the tracker closest to the calorimeter. Comparisons of low- 
background and high-background on-orbit data sets as well as comparisons of 7 -ray and CR- 
background Monte Carlo simulations have shown that these events suffer a higher background 
contamination fraction. 

Deposited charge veto: 7 rays convert into an electron-positron pair while most of the back- 
ground events involve a single charged particle. Therefore, we require the charge deposited 
in the first tracker layer following the interaction vertex to be > 1.5 times the value expected 
for a minimum ionizing particle, which is typically indicative that two particles crossed the 
layer rather than one by itself. 

Incidence angle veto: Events arriving from directions > 72° off the LAT boresight are rejected 
because there is increased CR background leakage for such highly inclined events. 


The new event class for the low-energy sample is denoted as P7REP_IGRB_L0 in the remainder 
of this manuscript to distinguish it from the publicly available standard event classes. The sky- 
averaged exposure of the P7REP_IGRB_L0 selection is 66 % of the exposure of the corresponding 
P7ULTRACLEAN selection for survey mode observations (see Figure [I]), when compared at the energy 
of maximum exposure. The estimated residual CR background rate is reduced by a factor of ~ 3 
around 200 MeV, where the background rate is highest. 


As a final step to define the low-energy sample, events from measured directions > 90° off 


the Earth zenith are rejected to limit contamination by photons from the Earth limb (Abdo et al. 
2009a). 


For the high-energy sample, we use a relaxed event selection compared to P7REP_IGRB_L0. 
At energies above 13 GeV, CR primaries in the form of protons and heavier nuclei dominate the 
background flux. The rejection power for CR nuclei is sufficient for this analysis if one requires 
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only the condition described above as ‘Tracker veto I’ in addition to the standard P7ULTRACLEAN 
event class definitions. We implement the ‘Incidence angle veto’ as for the low-energy event class. 

The standard P7ULTRACLEAN event classification rejects events for which the positions of the 
primary interaction vertex and the reconstructed shower maximum are separated by > 12 radiation 
lengths as measured along the shower axis. This selection criterion was introduced to reject CR 
events with bad shower reconstructions that would sometimes result in large apparent depths for 
the shower maxima — a strategy that works well for energies < 500 GeV, but removes a significant 
fraction of 7 rays > 500 GeV. Therefore this selection criterion is removed in the event selection for 
the high-energy sample. The very moderate increase in residual CR background arising from this 
removal is overcompensated by the ‘Tracker veto I’ condition that was introduced for this event 
class. 

The distinct classification scheme for the high-energy sample is denoted as P7REP_IGRB_HI in 
contrast to the P7REP_IGRB_L0 classification used for the low-energy sample. At high energies, the 
reconstructed arrival directions of CR-induced atmospheric 7 rays are confined to angles very close 
to the Earth limb (113° from the Earth zenith) due to the reduced width of the PSF (~ 0.1° above 
10 GeV). Therefore the zenith angle veto condition described above for the low-energy sample is 
modified to reject only photons from directions > 105° off the Earth zenith. 

The P7REP_IGRB_HI selection has a peak exposure of about 85% of the peak exposure of 
the standard P7ULTRACLEAN selection, and surpasses the P7ULTRACLEAN selection in acceptance 
> 700 GeV. 

Both new event selections were cross-validated against the standard P7ULTRACLEAN event 
selection by performing the analysis described below also on the P7ULTRACLEAN dataset. We 
obtain consistent results in the energy range in which we can perform the analysis even in the 
presence of the higher CR background of the P7ULTRACLEAN selection (400 MeV - 100 GeV). 
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Table 1. Event selection criteria for low-energy and high-energy samples, including modifications 

with respect to P7ULTRACLEAN 



Low-energy 

High-energy 

Data sample / event selection 

P7REP_IGRB_L0 

P7REP_IGRB_HI 

Add tracker veto I 

Y 

Y 

Add tracker veto II 

Y 

N 

Add deposited charge veto 

Y 

N 

Remove calorimeter shower maximum veto 

N 

Y 

Incidence angle veto 

> 72° 

> 72° 

Zenith angle veto 

> 90° 

> 105° 


Note. — See Section [Tl] for detailed descriptions of the event selection criteria. 
The low-energy and high-energy event samples are used to derive the IGRB inten- 
sity in the 100 MeV-13 GeV and 13 GeV-820 GeV energy ranges, respectively. 
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2.2. Instrument response functions 


New sets of dedicated instrument response functions (IRFs) were created for the low-energy 
(P7REP_IGRB_L0) and the high-energy (P7REP_IGRB_HI) event classes via Monte Carlo simulation of 
7 rays. The energy range of the new IRFs is 17.8 MeV to 1.78 TeV. 

Figure [l] shows the sky-averaged exposures obtained for the low-energy and the high-energy 
samples using the corresponding P7REP_IGRB_L0 and the P7REP_IGRB_HI IRFs, respectively. The ex- 
posure that would be obtained for the same observation period, but using the standard P7 JJLTRACLEAN 
event sample with IRFs P7REPJJLTRACLEAN_V15 is plotted for comparison. For the P7REP_IGRB_HI 
selection, there is an overall drop in exposure due to using part of the silicon tracker to veto 
charged particles. However, at the highest energies, especially > 300 GeV, this loss is increasingly 
counteracted by the removed shower maximum constraint in the P7REP_IGRB_HI class (see Section 


2.1). The P7REP_IGRB_L0 selection has a significantly lower average exposure than P7REP_IGRB_HI. 


This loss of exposure is acceptable at low energies where this event class is used since the IGRB 
measurement is not limited by statistics below tens of GeV. 

In-flight PSF corrections available for the IRFs corresponding to standard event classes have 
not been applied to the P7REP_IGRB_L0 and P7REP_IGRB_HI IRFs. The corrections were motivated 
by small differences observed in the PSF of the original (P7) on-orbit and simulated LAT data at 


energies > 1 GeV (Ackermann et al. 2012 f, 2013a). We verified directly that such small corrections, 


mitigated in the reprocessed data (Bregeon et al. 2013), do not significantly affect this analysis. 


This is expected since it is performed on a spatial grid of about 0.9 deg, considerably larger than 
the typical high-energy PSF. 


2.3. Residual cosmic-ray background 


Charged and neutral CRs misclassified as 7 rays by the multivariate event classification al- 
gorithms mimic an isotropic flux that is indistinguishable from the IGRB. In addition, genuine 7 
rays from the Earth’s atmosphere that have directional reconstruction errors sufficient to bypass 
the zenith angle veto become a source of apparent extraterrestrial emission over the full sky. In 
this work, the term ‘CR background’ includes CR-induced 7 rays from the atmosphere. 


Our estimation of the residual CR background event rate is based on Monte Carlo simulations of 
the relevant particle species in the near-Earth environment, namely CR nuclei and electrons, as well 
as their atmospheric secondaries. We simulate both CR backgrounds and signal 7 rays and extract 
the distributions for reconstructed event properties with the greatest discrimination power at low 
and high energies respectively: the multivariate event classifier output and the transverse shower 
size. The distributions for simulated background and signal are compared to the distributions for 
the flight data to quantify the level and associated uncertainty of the CR background. A detailed 
description of this method can be found in Ackermann et al. ( 2012 f). 
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To account for atmospheric 7 rays surviving the zenith angle veto, an updated phenomeno- 
logical model for the Earth emission based on LAT observations is included in the Monte Carlo 
simulation. Atmospheric 7 rays can bypass the zenith veto either by being reconstructed in the 
extreme tail of the PSF, or by entering from the back side of the instrument and being recon- 
structed as though coming from the front. Although such catastrophic mis-reconstructions are 
rare, the Earth emission is sufficiently bright that the expected event rate is non-negligible at en- 
ergies < 1 GeV (Bechtol 2012). For the stringent zenith angle selections used in this work, the 
residual contamination of atmospheric 7 rays is expected to be comprised primarily of back-entering 
events. The reconstructed directional distribution of back-entering atmospheric 7 rays in particular 
is approximately isotropic. 


Figure [2] shows the residual CR background rates as a function of reconstructed energy for the 
P7REP_IGRB_L0 and the P7REP_IGRB_HI classes. Note that the event energy is reconstructed under 
the hypothesis of a front-entering 7 ray and in general does not represent the actual energy for 
hadrons. At high energies, primary protons and electrons both contribute significantly to the CR 
background. Although protons are far more abundant than electrons in the environment of the 
LAT, there is also greater rejection power against protons since analysis of the shower shape in the 
calorimeter can be used to tag and remove protons in addition to the veto power obtained from the 
ACD. All contributions shown in Figure[2]have been adjusted from the raw Monte Carlo predictions 
based on event property comparisons between the simulated and flight data, as described above. 
The CR background contamination after re-scaling agrees with the raw predictions from Monte 
Carlo simulation to within 35%, depending on energy, and this maximum discrepancy is used as a 
measure of the systematic uncertainty in the CR background contamination. 


The full uncertainty in residual CR background rates, shown in Figure [2j has been calculated 
by adding systematic and statistical uncertainties in quadrature. For the P7REP_IGRB_HI event 
class at energies above 10 GeV (relevant for the high-energy sample), the statistical uncertainties 
are large due to the limited size of the simulated residual background sampk£j] Therefore, instead 
of using bin-by-bin estimates for the CR background rates in the high-energy IGRB analysis, the 
rates are obtained from a fit to the (re-scaled) simulated background rates between 10 GeV and 
820 GeV using a simple broken power law function in energy with a break at 50 GeV. A spline 
interpolation of the background rates is used as the CR background model below 10 GeV. 


4 The existing background rate estimates were derived from several million CPU hours of CR simulation. Significant 
gains in precision might be achieved in the future when more computing power becomes available. 
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Fig. 1. — Comparison of the sky-averaged exposure for the P7REP_IGRB_L0, P7REP_IGRB_HI, and 
P7ULTRACLEAN event selections. Thick lines indicate the respective energy ranges for which the 
P7REP_IGRB_L0 and P7REP_IGRB_HI event classes are used in this analysis. 
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Fig. 2. — Comparison of the rates of residual CR background events encountered in the 

P7REP_IGRB_L0 (left) and P7REP_IGRB_HI (right) event selections. The individual contributions 
from primary protons, electrons, secondary CR, and 7 rays produced in the atmosphere are shown 
based on the respective Monte Carlo predictions. The white crosses and gray boxes denote the 
total CR background contamination level including (gray boxes) and not including (white crosses) 
the systematic uncertainties. A band consisting of three black lines displays the model which is 
used for the level (thick line) and the uncertainty (thin lines) of the CR background in the IGRB 
analysis. 
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3. Derivation of the IGRB spectrum 


The spectrum of the IGRB is derived in a two-step procedure. First, the spectrum of the 
isotropic component is determined as part of a multi-component maximum likelihood fit of template 
maps to the observed LAT counts using the tools and method described in Ackermann et al. (2009). 
This isotropic component is attributed to the sum of the IGRB and misclassihed CR backgrounds 
(Section |2.3[ ). Second, the CR background contribution is subtracted from the isotropic component 
to obtain the IGRB spectrum. 


In the multi-component maximum likelihood fit, we create template maps containing the num- 
ber of LAT counts expected during the observation period for various contributions to the 7 -ray 
sky. Each template map is based on a model or measurements of the respective contribution. The 
7 -ray emission observed by the LAT is modeled using five template maps in addition to the isotropic 
emission and 403 point sources that are fitted individually. Two template maps describe the DGE. 
One map is used for the solar disk and inverse Compton (IC) emission associated with the solar 
radiation field. One map describes the local diffuse emission from Loop I and the Local Loop. A 
last map is used to describe contributions from established point sources that are not individually 
fitted. 


The normalization of each template is fitted individually for each energy bin in the energy 
range between 100 MeV and 13 GeV using the low-energy event sample. This fit is performed by 
maximizing the likelihood to obtain the observed number of counts in each pixel given our model. 
We denote this fit result as the ‘low-energy fit’. Above 13 GeV, the normalizations of the Galactic 
foreground templates are kept fixed over the full high-energy range, i.e. we use the spectral shape 
that is provided in the templates to model the foreground in addition to the spatial information. 
Only the normalizations of the isotropic template and the point sources are fitted for each individual 
energy bin above 13 GeV. To determine the fixed template normalization factors, we first fit the 
normalization of each Galactic foreground template in the six energy bins between 6.4 GeV and 
51 GeV using the same procedure as for the low-energy fit (the number of events above 51 GeV 
is too low to robustly fit all foregrounds individually in each energy bin). For each respective 
foreground template, the average normalization factor from these six energy bins is applied to all 
of the energy bins between 13 GeV and 820 GeV. The bin-by-bin fitting and average-value scaling 
procedures yield consistent spectral forms in the 6.4 GeV to 51 GeV range (cf. Figures 14 16, and 
[T 8 | Appendix |A|) , providing confidence to the extrapolation above 51 GeV based on the spectral 
model shape of each foreground component. Moreover, the shapes of the high-energy spectra of 
the dominant Galactic foreground contributions (i.e., the interactions of CRs with ISG and ISRF) 
can be derived quite robustly based on the well measured local CR electron and nucleon spectra. 
Modeling uncertainties are therefore small. 


After the Galactic foreground template normalizations are fixed to these average values, a 
multi-component maximum likelihood fit is performed using the high-energy event sample in the 
energy range between 13 GeV and 820 GeV, which we denote as the ‘high-energy fit’. 
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The template maps are binned on a HEALPi^grid of order 6 (~ 0.9° pixel size; Gorski et al 


2005) in Galactic longitude and latitude, and in 26 energy bins between 100 MeV and 820 GeV. 


Galactic diffuse emission dominates the intensity of the 7 -ray sky at low Galactic latitudes. This 
emission originates from the interactions of CRs with ISG and the ISRF. In the former case, 7 
rays are produced through the generation and decay of neutral pions, or through non-thermal 
bremsstrahlung; in the latter case, 7 rays are produced by the IC process. We consider the 7 - 
ray emission due to interactions with ISG and the ISRF separately in this analysis. The spatial 
distribution of 7 rays from interactions with ISG is well correlated to the column density of ISG 
along a given line of sight, whereas 7 rays created through interactions with the ISRF form a 
comparatively smooth emission component peaked at the Galactic center. We use the GALPROfj^] 
code ( Strong et al.pOOO Vladimirov et ah||2011 ) to obtain templates for these two components. A 
detailed discussion of the modeling of the Galactic diffuse emission will be presented in Section [4] 
and Appendix [A] We refer to the two template maps as the ‘Hi + Hil’ template and the ‘Inverse 
Compton’ template below. 

The second prominent contribution to the 7 -ray sky is from the individually resolved LAT 


7 -ray sources. We include the 403 sources from the second LAT source catalog (2FGL; Nolan et al. 


2012 ) from above and below the Galactic plane (| 6 | > 2 °) that are detected with a test statistic (TS) 


larger than 200 in that catalog as individual templates. The 1215 sources with TS values less than 
200 are added to a common template using the spectral information found in the 2FGL catalog. 


Additionally, we add a template for a source at the position of CRATES J231012— 051421 (Healey 


et al. 2007) after a localized excess in the residual map was found in a first iteration of this analysis 
at a position consistent with the CRATES source. Due to the difference in the observation time 
between this work and the 2FGL catalog (50 months vs. 24 months) and the intrinsic variability 
of many high-latitude LAT sources, both extra sources and changes in their time averaged spectra 
can be expected. However, no systematic search for additional sources too faint to be identified 
on the residual map was performed on the data samples used in this analysis. We note that due 
to the considerably more stringent event selection used in this work in comparison to the 2FGL 
catalog, the effective gain in exposure is much smaller than what is indicated by the difference in 
observation length. The threshold in 7 -ray flux for detecting sources in this sample would only be 
marginally lower than the flux threshold for the source sample in the 2FGL catalog. 

A third contribution to the 7 -ray sky that is included as a component in the likelihood fit 
is the 7 -ray emission related to the Sun. 7 rays are produced by CR collisions with the outer 
atmosphere of the Sun and by IC scattering of CR electrons off the solar radiation field. The solar 
disk and extended IC emission has been measured using the LAT (Abdo et ah] 2011). We use the 
Solar System Tools (Johannesson &; Orlando 2013) from the LAT ScienceTools (v9r30p0) to create 
a template for the time-averaged solar emission in our observation period which is based on this 


"http : / /healpix . jpl .nasa.gov/healpixSoftwareGetHealpix . shtml 
( http : / /galprop . stanf ord. edu 


-18- 


measurement. This template is denoted as ‘Solar Disk + IC’ throughout this work. To avoid bias 
at high energies where the solar spectrum is not well known by measurements, we do not use the 
±3° region around the ecliptic plane in the fit for energies above 13 GeV (in the energy range where 
the high-energy sample is used). The 7 -ray emission from the Moon (Abdo et al. 2012) has been 
neglected in this work since the Moon does not feature an extended IC contribution that could bias 
a measurement of the IGRB. 


Structures are seen in the diffuse 7 -ray emission that are correlated to Loop I (seen most 
prominently in the region of the North Polar Spur; Casandjian & Grenier [2009). Loop I is also 
bright in the 408 MHz radio continuum survey of Haslam et al. ( 1982) indicating a local overdensity 
of high-energy electrons or stronger magnetic fields in that region. A detailed investigation of the 
spectrum and spatial distribution of the 7 -ray emission from this region has not yet been performed, 
but is outside the scope of this paper. We therefore use a simple geometrical model (Wolleben 2007) 
for the synchrotron emission from Loop I and the Local Loop to generate a template for the 7 -ray 
emission from these structures (referred to as the ‘Loop I / Local Loop’ template). 


Systematic uncertainties associated with the foreground templates mentioned above and other 
foreground modeling choices, e.g., the optional inclusion of an additional template for the Fermi 
Bubbles, are discussed in Section |5.2| 

Certain regions in the vicinity of the Galactic plane have been masked and the corresponding 
pixels have not been used in the likelihood fit. The shape of the mask has been chosen to reduce 
systematic uncertainties connected to the Galactic diffuse foreground emission by excluding regions 
in which the column density of the ISG is not dominated by the atomic and ionized hydrogen in 
the vicinity of our solar system. Details of the definition of the mask are listed in Appendix [C| 
Figure [3] shows the integrated LAT counts above 100 MeV that are used for this analysis and the 
excluded regions. 
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Fig. 3. — Map of counts observed by the Fermi LAT above 100 MeV using a Mollweide projection 
in Galactic coordinates with a pixel scale of ~ 0.9°. The color scale is logarithmic. Overlaid is the 
mask used in this analysis to exclude regions from the template fitting procedure (see Appendix [C] 
for details). 
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4. Foreground diffuse Galactic emission models 


At high Galactic latitudes (\b\ > 10°) the ISG is dominated by atomic gas clouds within a 
few hundred parsecs, a range in which we do not expect a significant change in the density of 
the CRs that interact with these clouds (Ackermann et al. 2013b). Therefore the spectrum and 
spatial distribution of the 7 rays arising from the CR interactions with atomic gas are relatively 
well constrained by measurements of the gas cloud distributions and direct measurements of local 
CR spectra. The proximity of ISG seen at high Galactic latitudes permits the use of a single model 
template as opposed to multiple templates describing the ISG at various distances, as is typical 
for other DGE studies. Only a small fraction of the ISG-related 7 -ray emission arises from ionized 
hydrogen gas, which has a larger scale height than the atomic gas and a less well known distribution 
on the sky. 


The spectrum and intensity distribution of the IC-related 7 -ray emission on the sky can only 
be predicted from a global modeling of CR propagation and interaction in the Galaxy. It is highly 
dependent on the injection and propagation of electrons in the Galactic plane and into the Galactic 
halo. It further depends on the spatially varying spectrum of the ISRF, which is more uncertain 
than the distribution of the ISG. 


4.1. Reference foreground models 


An extensive study comparing LAT data to GALPROP-based predictions of the DGE has 


been published in Ackermann et al. (2012b). However, that study was restricted to a selected set of 


propagation and CR injection scenarios. In particular, it was assumed that the diffusion coefficient 
is constant throughout the Galaxy and the source population that injects the electrons is the same 
as the source population that injects the nuclei. These choices are well motivated by Occam’s 
razor, but we demonstrate in Appendix |A| that template maps for the IC emission that are derived 
from DGE models of this type lead to inconsistencies when used in our multi-component likelihood 
fit. For example, the spectrum of the IC emission predicted by the model does not match the 
spectrum obtained in the fit to the LAT data. Such a mismatch is critical for the IGRB study if 
it originates from an inaccurate model of the IC intensity distribution on the sky. In this case, the 
isotropic template could partially compensate for the inaccuracies of the IC intensity distribution 
and thereby lead to a biased IGRB measurement. 

We therefore extend the study of foreground models to include two reference models for prop- 
agation and injection scenarios with more degrees of freedom than those considered in Ackermann 


et al. (2012b). One allows for distinct populations injecting CR electrons and nuclei. The other al- 


lows a variation of the diffusion coefficient with radial distance from the Galactic center and height 
above the Galactic plane. These two reference models are described in more detail in Appendix [A] 
We denote them as foreground models ‘B’ and ‘C’ respectively, to distinguish them from foreground 


model ‘A’ that is derived from the class of DGE models studied in Ackermann et al. (2012b). The 
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principal features of and differences between the three reference foreground models are summarized 
in Table [2j 


Our main concern in this work is to investigate whether the fitted spectral features of the 
IGRB depend on the specific type of foreground model chosen. It is not the aim or scope of this 
work to perform a quantitative study of whether one of the alternative foreground model classes 
matches the LAT data better than another. For simplicity, we use model A as a baseline for the 
purpose of quoting certain results and testing variations of DGE model parameters. However, we 
do not view model A as canonical or preferred over the other models. 


All three foreground models assume diffusive CR transport with re-acceleration in the interstel- 
lar medium. The diffusion coefficient has a power-law dependence on rigidity with index S = 0.33, 
as expected from a Kolmogorov spectrum of magnetic turbulences. CR propagation and injection 
parameters within each model are chosen to obtain good agreement between the predicted local 
spectra of interstellar CRs and actual CR measurements after solar modulation effects have been 
taken into account. In particular, we require good agreement with the measured proton and helium 
spectra, the electron spectrum, and the B/C and 10 Be/ 9 Be ratio^J 


The distribution of H 2 and Hi gas in the Galaxy is modeled based on the microwave survey of 
|Dame et al. (2001) and the radio survey of Kalberla et al. (2005), respectively (see also Appendix 
0. Details regarding the gas distribution modeling can be found in Ackermann et al.| (2012b|). An 
analytic model is used for the distribution of ionized hydrogen in the Galaxy (Gaensler et al. 2008). 
We use the ISRF model introduced in Porter et al. (2008), which is available within GALPROP. 
We take into account the anisotropy of the ISRF by calculating for 192 uniformly distributed lines- 
of-sight the ratio between the predicted IC emission from a full anisotropic calculation and the 


7 One notable exception is the electron spectrum for foreground model B that is tuned to reproduce the observed 
IC emission spectrum rather than the measured CR electron spectrum. A detailed description of the various injection 
and propagation parameters for all three foreground models can be found in Appendix [Al 


Table 2. Comparison of Benchmark Galactic Foreground Models 


Foreground 

Main features and differences with respect to other DGE models 

Model A 

Sources of CR nuclei and electrons trace pulsar distribution; 
constant CR diffusion coefficient and re-acceleration strength through Galaxy 

Model B 

Additional electron-only source population near Galactic center, 
these electrons are responsible for majority of IC emission; 
local source of soft CR electrons needed to explain CR electron spectrum at Earth below 20 GV 

Model C 

Sources of CR nuclei and electrons more centrally peaked than pulsar distribution; 

CR diffusion coefficient and re-acceleration strength vary with Galactocentric radius and height 


Note. — See Appendix [A] for a more detailed description of these three reference DGE models. 
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prediction assuming that the ISRF is isotropic. This set of ratios is then interpolated and applied 
as a multiplicative correction to all the generated IC templates. 


It has been shown that interstellar dust can trace gas that is not seen in Hi or CO surveys 
(Grenier et al. 12005). Therefore we use the E (B — V ) visual reddening map provided by Schlegel 


et al. (1998), a tracer of the interstellar dust column density, to estimate the total ISG column 


density along a line-of-sight. We use the procedure described in Ackermann et al. (2012b) to 
obtain a conversion factor between the magnitude of reddening in the E (B — V) map and the Hi 
gas column density (denoted as Hi-to-dust ratio below). The procedure uses a linear regression 
between the E (B — V) map and the Hi and CO surveys to obtain the Hi-to-dust ratio in regions 
of the sky where E (B — V)< 5 mag. The fit depends slightly on the spin temperature Tg that one 
assumes to correct for the opacity of the 21 cm-line in the Hi surveys. We find a Hi-to-dust ratio 
of 7.9 x 10 21 cm' 2 mag -1 assuming the widely used value of Tg= 125 K for the spin temperature 
(e.g., Kulkarni & Heiles 1988), and use this Hi-to-dust ratio in our analysis. Note that this value 
is different from the value used in Ackermann et al. (2012b) where the higher spin temperature of 
Tg = 150 K was used for deriving this conversion factor. The lower spin temperature used here 
leads to smaller residuals between the E (B — V ) and the Hi and CO survey derived gas column 
densities at high Galactic latitudes that are relevant for the IGRB analysis. 


4.2. Additional foreground models used for systematics studies 


In addition to our three reference foreground models, we consider further variations of fore- 
ground models to assess the systematic uncertainties for the derivation of the IGRB that are related 
to the modeling of the DGE. Specifically, we study variations of the size of the CR halo between 
4 kpc and 10 kpc, a variation of the distribution of CR sources in Galactocentric radius, a model 
where CRs are not re-accelerated in the Galaxy, models with an extra foreground template for the 
Fermi Bubbles, a higher radiation field in the Galactic Bulge, a lower random Galactic magnetic 
field than in the default models, and a variation of the Hi-to-dust ratio by 10%. 


Foreground model A serves as the baseline model for these variations. To assess the stability 
of the IGRB measurement with respect to assumptions regarding the halo size and the CR source 


distribution we use three models discussed in Ackermann et al. (2012b) chosen to cover the extreme 


values for the CR halo size (4 kpc vs. 10 kpc) and radial source distribution (traced by pulsars vs. 
traced by SNRs) in the range of models studied there. 


Strong et al. (2011) suggest that CR propagation models where CRs are not re-accelerated in 


the Galaxy (so called plain diffusion models) describe the synchrotron emission observed from the 
Galaxy better than present models with re-acceleration. We therefore include such a plain diffusion 
model in our investigations. 


Dobler et al. ( 2010| ) and Su et al. (2010) have noted the existence of large-scale structures of 
residual diffuse 7-ray emission above and below the Galactic center region that became well known 
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as the ‘ Fermi Bubbles’. Although well established as significant features, the Bubbles were not 
included in the reference foreground models. The Fermi Bubbles have been studied exclusively as 
a residual after subtracting a model of the DGE and isotropic diffuse emission. No template for 
their shape has yet been derived from independent observations. Using a template derived from 
a strictly empirical excess of 7 -ray emission might lead to a bias in the other components of the 
fit, including the isotropic template. Since neglecting the emission from the Fermi Bubbles might 
bias the fit as well, we tested the effects of including template maps for the Fermi Bubbles in 
the multi-component fit. Two models for the intensity distribution of the Fermi Bubbles on the 
sky were tested. The first is a simple geometric template used for the investigation of systematic 


uncertainties in studies of the 7 -ray emission from SNRs ( 

de Palma et al. 

2013 

). The second is a 

template derived from the residual 7 -ray emission ( 

Fermi-LAT Collaboration 

2014). 


Our knowledge of the ISRF in the inner Galaxy is limited. We therefore repeat the IGRB fit 
against an ISRF model that assumes a factor ten higher stellar luminosity in the Galactic Bulge. 
Such a model is still compatible with constraints derived from observations of the ISRF in the solar 
neighborhood. 


We also test the impact of the assumed random magnetic field strength by generating a fore- 
ground model with a lower random magnetic field strength of 3 /iG in the solar neighborhood, to 
compare with the reference models A and B which use a value of 7.5 /rG for the random magnetic 
field there (see in this context Appendix | A. 1[). 


The difference in the Hi-to-dust ratio between the spin temperature value of Tg = 150 K 
adopted in Ackermann et al. (2012b) and the widely used spin temperature value of Tg = 125 K 
that was used in this work is of order 10 %. We test foreground models with variations of the 
Hi-to-dust ratios by ±10 % to determine the impact of this parameter on the IGRB measurement. 


The impact of all described variations in modeling the foreground DGE on the spectrum of the 
IGRB is discussed in Section 5.2 together with a model-independent study to assess the impact of 


un-modeled residuals in the foreground emission. 
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5. Results 


5.1. IGRB spectrum 


The likelihood fitting technique introduced in Section [3] is used to derive the spectrum of 
the isotropic emission for the three different DGE foreground models described in Section |4j The 
residual particle background contamination is subtracted from the isotropic component to obtain 
the spectrum of the IGRB. 


Figures[4j[5j and[6]show the results of the fits using foreground models A, B, and C, respectively. 
Each figure displays the average high- latitude (|6| > 20°) intensities attributed to the isotropic 
emission (IGRB plus CR background), the individual sources, two DGE components, the solar 
emission, and the local foreground templates. The sum of these intensities is compared to the 
average 7-ray intensity observed by the LAT. A separate graph shows the contributions of the 
IGRB and of the residual CR contamination to the isotropic emission. Since the isotropic and 
IGRB intensities in the highest energy band are compatible with zero within the lcr uncertainty 
range we quote upper limits in that energy banc0 The error bars displayed for the individual 
components in the three figures include the statistical uncertainty and the systematic uncertainty 
of the effective area parametrization (Ackermann et al. 2012f) added in quadrature. The error bars 
for the IGRB component additionally contain the systematic uncertainty due to subtracting a not 
perfectly known CR background contamination, also added in quadrature. 


The IGRB intensities corresponding to foreground models A, B, and C are compared in Figure 
[Tj Numerical values for the IGRB intensities per energy band when using foreground model A 
are presented in Table [3} Intensities for other foreground models can be found in the electronic 
supplementary material to this article. The IGRB intensity shows a clear cutoff at high energies, 
independent of the foreground model. A x 2 regression of the IGRB spectrum using a power law 
with exponential cutoff (PLE) spectral model of the form 


IE = /lM (ioofiev) “ P (l£) (1) 

results in low x ' 2 values for all three foreground models and can therefore be considered suitable to 
characterize the IGRB spectrum. We further try to fit the IGRB intensities with a single power 
law (PL) and a smoothly connected broken power law (BPL). The fit parameters for the PLE 
model, as well as the x 2 values for all fitted spectral hypotheses are summarized in Table [4| The 
X 2 values for the BPL and the PLE spectral models are similar enough that the two hypotheses are 
indistinguishable in the energy range observed. We prefer to quote fitted parameter values for the 


8 Statistical uncertainties on the isotropic emission have been calculated using the MINOS algorithm of the MINUIT 
minimization package ( |James fe Roos|1975| ). The position of the upper limit corresponds to the upper bound of the 
la uncertainty interval. 
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PLE model given its lower number (three vs. four) of free parameters. The PL model is disfavored 
independent of the foreground model based on the high y 2 values of 88 (foreground model A), 
151 (foreground model B), and 106 (foreground model C) for 23 degrees of freedom. Note that the 
y 2 value cannot be easily interpreted in terms of a significance for the agreement between spectral 
model and data because the error bars of the IGRB spectrum are systematics dominated over most 
of the energy range, and therefore correlations between bins are expected. These correlations are 
also responsible for the rather small y 2 values (when compared to the number of degrees of freedom) 
for the PLE and BPL spectral models. 

The large difference in y 2 between the PL and the PLE models even when neglecting bin-to- 
bin correlations can still be interpreted as a robust evidence against a simple power-law spectrum. 
For the benchmark models this y 2 difference is larger than 61 at just one added degree of freedom 
in the model. We also calculated the y 2 differences between the PL and PLE models for the 
additional foreground models used in the investigation of the foreground related systematics that 
are summarized in Table 4. The y 2 difference between the PL and PLE models is 45 or larger in 
all of these foreground variations. 

Finally we calculated the y 2 difference between PL and PLE models if we add the foreground 
model related systematic uncertainties to the instrument related uncertainties. Since we do not 
know the correlations between the bins introduced by the systematic errors we adopt a worst-case 
scenario here assuming that the dominant fraction of the systematic error is fully correlated at low 
energies and anti-correlated between low and high energies (E >300 GeV). Such a hypothetical 
anti-correlation in the systematics might artificially enhance indications of a cutoff. Even in this 
worst-case scenario we still find y 2 differences exceeding 25 between the PL and PLE scenarios for 
our 3 benchmark models. 

Residual maps of the relative deviations in intensity between model and data in different regions 
of the sky can be found in Appendix [B] None of the models considered is a perfect description of 
the data and large-scale residuals at the 25% level appear in many parts of the sky. No Galactic 
foreground model is unambiguously preferred over the others based on the level and distribution of 
the residual 7 -ray emission. Prominent residuals include the Fermi Bubbles, but other un-modeled 
residuals appear in different parts of the sky. Also, there seems to be less intensity observed in the 
south polar region than its northern counterpart, a feature that cannot be modeled with the classes 
of foreground models considered here which exhibit a symmetrical CR density about the Galactic 
plane by construction. 

In the next section, we present a study of the general impact of such un-modeled residuals on 
the IGRB spectrum. The study is not restricted to the Fermi Bubbles, but applies to all features 
in the residual maps. 
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5.2. Systematic uncertainties from foreground modeling 

The fitting procedure is applied to the variants of the foreground models introduced in Section 


parameter and y 2 values obtained for the IGRB when using each of the foreground model variants. 
For most of the variant models, the general shape of the IGRB spectrum is not affected; a simple 
power law with exponential cutoff remains a valid parametrization. The single exception is the 
effect from changing the source distribution. When using the distribution of SNRs as a tracer of 
the CR source density, a second apparent spectral feature manifests as a dip in the IGRB spectrum 
at a few GeV. This dip in the IGRB compensates the higher intensity of the IC emission with 
respect to other foreground models. A good fit of the shape of the IGRB spectrum can be obtained 
only by describing it as a sum of two components. We choose two components each having a power- 
law plus exponential cutoff spectrum to describe the resulting IGRB spectrum for this foreground 
model version. 

A second investigation of systematic uncertainties is aimed at the residuals visible in the maps 
shown in Appendix [B] We study the variations of differences between LAT data and our model 
within four very high-latitude (\b\ > 60°) overlapping regions in which the IGRB is responsible for a 
substantial fraction of the total 7 -ray emission. We use the Galactic north pole, the Galactic south 
pole, the | 6 | > 60° region facing the inner Galaxy, and the | 6 | > 60° region facing the outer Galaxy. 
For each of these regions, we calculate the spectral residual and the renormalization factor for the 
IGRB needed to obtain the best overall agreement between model and data in the corresponding 
region. The resulting adjustment factors are between 0.7 and 0.95 depending on the region and 
foreground model. Note that all adjustment factors are < 1, i.e., the observed 7 -ray intensity at 
| 6 | > 60° is lower than that predicted by the models. This overestinration of the 7 -ray intensity 
could arise either from an overestinration of the Galactic foreground at high latitudes or from an 
overestimation of the isotropic intensity, and therefore represents a systematic uncertainty for our 
analysis. 

Including this second study, we find that the normalization of the IGRB intensity i>ioo varies 
by +15%/— 30% with respect to foreground model A depending on our assumptions about the DGE 
foreground, while the spectral index varies between 2.26 and 2.34, and the cutoff energy between 
206 and 374 GeV. Summarizing the results above, the systematic uncertainty in the IGRB spectrum 
associated with foreground modeling is shown in Figure [7} 


4.2 in the same way as for the benchmark models A, B, and C. Table [5] summarizes the spectral 


5.3. Total EGB intensity 

The total EGB as defined in Section [I] is a quantity independent of the instrument and obser- 
vation time. For the purpose of this analysis, we call the sum of the IGRB and the sky-averaged 
intensity of high-latitude \b\ > 20° resolved LAT sources the total EGB. We note that this definition 
of the total EGB formally includes a small fraction of Galactic sources, e.g., millisecond pulsars. 
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However, the contributions of known high-latitude pulsars extracted from the second LAT pulsar 
catalog (Abdo et al. 2013) is less than 5% of the total EGB anywhere in the measured energy 
range, well below the level of systematic uncertainty inherent to the IGRB measurement. We fur- 
ther check the variation of the total intensity when varying the latitude threshold from \b\ > 20° to 
|6| > 40° to probe a possible Galactic source contamination. We find that the derived intensities 
are consistent within the respective uncertainties for each latitude threshold. 


Figure [8] compares the total EGB derived for foreground models A, B, and C, respectively. 
Numerical values for the total EGB intensities per energy band are given in Table [3j Again, we 
use a x 2 -regression to test different functional parametrizations of the spectrum. The best-fit 
parameters for a fit with a power law with an exponential cutoff as well as y 2 values for all tested 
spectral models are summarized in Table [6j For the total EGB we find similarly as for the IGRB 
that a power law with an exponential cutoff describes the spectral shape significantly better than 
an unbroken power law. The cutoff energy is higher for the total EGB than for the IGRB. As in 
the case of the IGRB, we cannot distinguish an exponential cutoff spectral model from a broken 
power law. Results of the spectral fits are summarized in Table [6j 

Systematic uncertainties in the total EGB spectrum arising from modeling the Galactic fore- 
ground are indicated by the shaded band in Figure [8j constructed using the identical methods 
described in Section lT2l for the IGRB. 
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Fig. 4. — Results of the IGRB fit for foreground model A. Average intensities for Galactic lati- 
tudes \b\ > 20° are shown. Left: Intensities attributed to the different foreground templates, the 
isotropic emission, and the individually resolved sources in the multi-component likelihood fit. The 
isotropic emission and the individually resolved sources are fitted separately in each energy bin. All 
other components are fitted individually in each energy bin below 13 GeV and included with fixed 
normalizations above this energy. The total intensity obtained from the IGRB fit is compared to 
the total intensity observed by the LAT. Error bars include statistical errors as well as systematic 
errors from the uncertainty in the LAT effective area parametrization. See text for details. Right: 
IGRB and CR background contributions to the isotropic emission. The line indicates the best-fit 
IGRB spectrum with a power-law plus exponential cutoff spectral model. Spectral parameters are 
given in the legend. The total intensities obtained from the fit and measured by the LAT are shown 
with the CR background contributions subtracted in this graph. 
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Fig. 5. — Results of the IGRB fit for foreground model B. See Figure [4] for legend. 




Fig. 6. — Results of the IGRB fit for foreground model C. See Figure [4] for legend. 
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Fig. 7. — Comparison of the derived IGRB intensities for different foreground (FG) models. The 
error bars include the statistical uncertainty and systematic uncertainties from the effective area 
parametrization, as well as the CR background subtraction (statistical and systematic uncertainties 
have been added in quadrature). The shaded band indicates the systematic uncertainty arising from 
uncertainties in the Galactic foreground: the IGRB intensity range spanned by the three benchmark 
models, the variants described in Section [4.2| and the normalization uncertainties derived from the 


high-latitude data/model comparison. See Section 5.2 for details. 
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Fig. 8. — Comparison of the total EGB intensities for different foreground models. The total EGB 
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Table 3. IGRB and total EGB intensities 


Energy range 

IGRB 

FG model uncert. 

Total EGB 

FG model uncert. 

Sources | 6 | > 20° 

[GeV] 

intensity 

on IGRB 

intensity 

on total EGB 



[cm - 2 s — 1 sr — 1 ] 

[cm - 2 s — 1 sr — 1 ] 

[cm - 2 s — 1 sr — 1 ] 

[cm 2 : 

s^sr- 1 ] 

[cm - 2 s — 1 sr — 1 ] 

0.10 - 0.14 

( 2.8 ± 0 . 6 ) x 10“ 6 

+ 01 y ID -6 
- 0.9 x iU 

(3.7 ± 0.6) x 10 ~ 6 

+ 0.3 

- 1.2 

X 10 ~ 6 

(9.0 ± 1.6) x 10 “ 7 

0.14 - 0.20 

(1.7 ±0.4) x 10 - 6 

tg;£ x 10-6 

(2.3 ±0.4) x IO " 6 

+ 0.1 

- 0.8 

X 10 -6 

( 6.2 ± 1 . 0 ) x 10 -7 

0.20 - 0.28 

(1.1 ± 0.3) x 10 ~ 6 

±S:i x w -6 

(1.5 ±0.3) x 10 ~ 6 

+ 0.1 

- 0.5 

x 10 -6 

(4.1 ± 0.6) x IO " 7 

0.28 - 0.40 

(6.7 ± 2.0) x 10 -7 

t°A X 10 -7 

(9.7 ± 2.0) x 10 ~ 7 

+ 1.2 
- 3.2 

x 10 -7 

(3.0 ±0.4) x 10 -7 

0.40 - 0.57 

(4.5 ± 1.0) x 10 “ 7 

t?(e X 10 - 7 

(6.7 ± 1.0) x IO ” 7 

+ 0.9 

- 2.2 

x 10“ 7 

( 2.2 ± 0 . 2 ) x 10“ 7 

0.57 - 0.80 

(3.3 ±0.4) x 10 “ 7 

X 10_7 

(4.9 ±0.4) x IO ’ 7 

+ 0.7 

- 1.6 

x 10 -7 

( 1.6 ± 0 . 1 ) x IO " 7 

0.80 - 1.1 

(1.9 ±0.2) x 10 “ 7 

tg;* x 10“ 7 

(3.0 ±0.2) x IO ’ 7 

+ 0.5 

- 1.0 

x 10 -7 

( 1.1 ± 0 . 1 ) x IO " 7 

1.1 - 1.6 

( 1.1 ± 0 . 1 ) x 10 -7 

tS:I x io- 7 

(1.8 ±0.1) x IO ’ 7 

+ 0.4 

- 0.6 

x 10“ 7 

(7.1 ± 0.7) x IO " 8 

1.6 - 2.3 

( 6.0 ± 0 . 8 ) x 10~ 8 

±2.7 x 10 “ 8 

( 1.1 ± 0 . 1 ) x 10~ 7 

±0.3 

x 10 -7 

(4.8 ±0.5) x IO " 8 

2.3 - 3.2 

(3.9 ±0.4) x 10 “ 8 

+1 x 10“ 8 

(6.9 ±0.5) x 10 ~ 8 

+ 1.9 
- 2.1 

x 10“ 8 

(3.0 ±0.3) x IO " 8 

3.2 - 4.5 

(2.3 ±0.3) x 10 ~ 8 

+? x 10“ 8 

(4.2 ±0.4) x 10 ~ 8 

+ 1.3 
- 1.2 

x 10“ 8 

(1.9 ±0.2) x IO " 8 

4.5 - 6.4 

(1.5 ± 0.2) x 10 ~ 8 

±6 x ! 0~ 8 

(2.6 ±0.3) x 10 ~ 8 

+ 0.8 

- 0.7 

x 10 -8 

( 1.1 ± 0 . 1 ) x 10“ 8 

6.4 - 9.1 

(9.6 ± 1.5) x IO " 9 

tit x io- 9 

(1.7 ± 0.2) x 10 ~ 8 

±0.5 

x 10“ 8 

(7.3 ±0.9) x IO " 9 

9.1 - 13 

(7.6 ± 1.0) x 10 “ 9 

t+ x IO " 9 

( 1.2 ± 0 . 1 ) x 10~ 8 

±0.3 

x 10 -8 

(4.4 ± 0.5) x IO " 9 

13 - 18 

(4.0 ±0.5) x 10 “ 9 

+f° x IO -9 

( 6.8 ± 0 . 6 ) x 10~ 9 

+ 2.0 

- 1.9 

x 10“ 9 

(2.7 ± 0.3) x IO " 9 

18 - 26 

(2.6 ±0.3) x IO " 9 

- 0.7 X HT 9 

(4.4 ±0.4) x 10 -9 

± 1.2 

x 10 -9 

(1.7 ± 0.2) x IO " 9 

26 - 36 

( 1.6 ± 0 . 2 ) x 10 -9 

to I X 10 - 9 

(2.7 ± 0.2) x 10 ~ 9 

±0.7 

x IO " 9 

( 1.1 ± 0 . 1 ) x 10 -9 

36 - 51 

( 1.1 ± 0 . 1 ) x 10~ 9 

to.3 X 10 ~ 9 

( 1.8 ± 0 . 2 ) x 10~ 9 

+ 0.4 

- 0.5 

x 10“ 9 

(7.3 ± 0.9) x IO ” 10 

51 - 72 

(6.3 ±0.8) x 10 - 10 

t, ° X IO - 10 

( 1.1 ± 0 . 1 ) x 10~ 9 

+ 0.2 

- 0.3 

x 10“ 9 

(4.5 ± 0.6) x IO ” 10 

72 - 100 

(3.6 ±0.5) x IO ” 10 

+J X IO - 10 

(6.2 ±0.6) x IO ” 10 

+ 1.1 
- 1.7 

x IO -10 

(2.6 ± 0.3) x IO ” 10 

100 - 140 

(1.5 ±0.3) x 10 ~ 10 

to 4 X 10 - 10 

(3.1 ± 0.4) x IO ” 10 

+ 0.5 

- 0.9 

x IO ” 10 

(1.5 ± 0.2) x IO ” 10 

140 - 200 

(9.8 ±2.0) x 10“ u 

tn x io- n 

(1.9 ± 0.3) x IO ” 10 

+ 0.3 

- 0.5 

x 10“ 10 

(9.3 ± 1.6) x IO ” 11 

200 - 290 

( 4 . 7 + 3 ) x 10~ n 

t, 2 X IO” 11 

(8.9 ± 1.7) x IO” 11 

+ 1.3 
- 2.4 

x IO" 11 

(4.1 ± 1.0) x IO” 11 

290 - 410 

(3.2 T i ( q ) x 10~ n 

tg; 1 x io- 11 

(6.3 ± 1.3) x 10” 11 

+ 0.9 

- 1.7 

x IO -11 

(3.0 ± 0.8) x IO” 11 

410 - 580 

(7.3 ±|;X ) x 10~ 12 

t 8 ;| x io -12 

(2.1 1 o(s ) X 10- 11 

+ 0.4 

- 0.5 

x IO -11 

(1.3 ± 0.6) x IO” 11 

580 - 820 

< 2.3 x 10” 12 


(9.7 ± 6.0) x IO” 12 

+ 2.3 

- 2.8 

x IO -12 

(9.0 ± 5.2) x 10“ 12 


Note. — Measured intensities of the IGRB, the total EGB, and the identified sources (|6| > 20°) per energy band, when 
using model A to describe the Galactic foreground. Uncertainties arising from foreground (FG) modeling are given in separate 
columns. Digitized versions of this table and the corresponding results for foreground models B and C are available in the online 
supplementary materials. 



Table 4. 


Results of the parametric fit of the IGRB 


Foreground 

model 

hoo 

[MeV - 1 cm - 2 s - 1 sr -1 ] 

7 

-Ecut 

[GeV] 

-f >100 

[cm -2 s -1 sr _1 ] 

X 2 /ndof 
(PLE a ) 

X 2 /ndof 
(PL a ) 

X 2 / ndof 
(BPL a ) 

Model A 

(0.95 ± 0.08) x 10 “ 7 

2.32 ±0.02 

279 ± 52 

(7.2 ±0.6) x 10 ~ 6 

13.9/23 

87.5/24 

13.5/22 

Model B 

(1.12 ± 0.08) X 10 “ 7 

2.28 ± 0.02 

206 ± 31 

(8.7 ± 0.6) x lO " 6 

7.9/23 

151. /24 

10 . 6/22 

Model C 

(0.78 ±0.07) x 10 “ 7 

2.26 ± 0.02 

233 ± 41 

( 6.2 ± 0 . 6 ) x 10" 6 

10.7/23 

106.5/24 

11.3/22 


a PLE = power-law plus exponential cutoff ; PL = power-law ; BPL = broken power-law 


Note. — Parameters obtained from a parametric fit of the IGRB spectrum. Intensity iioo? spectral index 7 and cutoff 
energy E cu t for a fit of the observed spectrum with the function given in Equation [l] are shown in columns 2-4. The integrated 
IGRB intensity above 100 MeV, />ioo 5 is found in column 5. A comparison of the x 2 / n dof values between the fit with the 
function in Equation [l] and alternative spectral models is given in columns 6 - 8 . The x 2 values include systematic uncertainty. 
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Table 5. Impact of foreground model variations on IGRB spectral parameters 


Model variation 

boo 

7 

E C ut 

X 2 /ndof 


[MeV _1 cm _2 s _1 sr _1 ] 


[GeV] 



10 x ISRF in Galactic Bulge 
Random magnetic field strength 3 
Plain diffusion model 


Model 38 a (Pulsars trace CR sources, 4kpc CR halo) 
Model 62 a (Pulsars trace CR sources, lOkpc CR halo) 
Model 06 a ’ b (SNR trace CR sources, 4kpc CR halo) 
Fermi Bubbles template A 
Fermi Bubbles template B 
Hi-to-dust ratio +10% 

Hi-to-dust ratio —10% 


(0.96 ±0.08) X 

10- 7 

2.31 ± 0.02 

273 ± 50 

13.4/23 

(0.93 ±0.09) x 

10- 7 

2.33 ±0.03 

257 ± 56 

14.1/23 

(0.91 ± 0.08) x 

10- 7 

2.31 ± 0.02 

264 ± 52 

13.4/23 

(1.06 ±0.09) x 

10- 7 

2.33 ±0.02 

367 ± 75 

15.0/23 

(0.89 ±0.08) x 

10- 7 

2.31 ± 0.02 

374 ± 77 

16.7/23 

(0.56 ±0.07) x 

10- 7 

2.27 ± 0.03 

399 ± 92 

55/23 

(1.02 ± 0.09) x 

10- 7 

2.32 ±0.02 

229 ± 44 

13.0/23 

(1.05 ± 0.09) x 

10- 7 

2.31 ± 0.02 

244 ± 42 

13.7/23 

(1.09 ±0.09) x 

10- 7 

2.34 ± 0.02 

280 ± 52 

12.0/23 

(0.82 ± 0.08) x 

IQ" 7 

2.30 ± 0.03 

274 ± 51 

17.5/23 


a from |Ackermann et alT] ( |2012b[) 

b Model is not well fit by a simple power-law with exponential cutoff spectral hypothesis. A two-component fit is used in- 
stead for this model for the evaluation of the foreground model systematics. The parameters obtained from this fit are I^ Q = 
0.69 x 10 -7 GeV _1 cm _2 s _1 sr _1 , 7^°) =1.74, E^ t = 0.60 GeV for the first, and i[q q = 0.16 x 10 -7 GeV _1 cm _2 s _1 sr _1 , 
7^) = 2.02, = 183 GeV for the second component. 


Note. — Parameters obtained from fits of the IGRB spectrum for variants of the DGE foreground model. Specific 
intensity /ioo 5 spectral index 7, and cutoff energy E cu t for a fit of the observed spectrum with the function given in 
Equation are shown in columns 2-4. The x 2 / n dof values of the fit are shown in column 5. 


Table 6. Results of the parametric fit of the total EGB 


Foreground 

model 

boo 

[MeV~ 1 cm _2 s _1 sr _1 ] 

7 

E C ut 

[GeV] 

-G100 

[cm -2 s — 1 sr — 1 ] 

X 2 / ndof 
(PLE a ) 

X 2 / ndof 
(PL a ) 

X 2 /ndof 
(BPL a ) 

Model A 

(1.48 ±0.09) X lO" 7 

2.31 ± 0.02 

362 ± 64 

(1.13 ±0.07) x 10~ 5 

11.0/23 

72.4/24 

10.5/22 

Model B 

(1.66 ±0.09) x 10" 7 

2.28 ±0.01 

267 ± 37 

(1.29 ±0.07) x 10~ 5 

13.5/23 

130./24 

11.3/22 

Model C 

(1.28 ±0.08) x 10" 7 

2.30 ±0.02 

366 ± 71 

(0.98 ±0.06) x 10~ 5 

6.9/23 

91.1/24 

7.7/22 


a PLE = power-law plus exponential cutoff ; PL = power-law ; BPL = broken power-law 


Note. — Parameters obtained from fits of the total EGB. Intensity /ioo 5 spectral index 7 and cutoff energy E cu t for a fit of 
the observed spectrum with the PLE function given in Equation [I] are shown in columns 2-4. The integrated IGRB intensity 
above 100 MeV, />ioo 5 is found in column 5. A comparison of the x 2 / n d°f values between the fit with the function in Equation 
^and alternative spectral models is given in columns 6-8. The % 2 values include systematic uncertainties. 
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6. Discussion and conclusions 


We have refined the measurement of the LAT IGRB intensity relative to the analysis of Abdo 


et al. (2010b), which was based on 10 months of LAT observations, now using 50 months of 


accumulated data. The measurement lower bound has been extended from 200 MeV to 100 MeV, 
and we report the first IGRB measurement with any instrument between 102 GeV and 820 GeV. 

The updated LAT IGRB spectrum remains consistent with a featureless power law between 
100 MeV and 100 GeV, and there is now strong evidence for a high-energy cutoff feature. The 
spectrum is well described by a power law with an exponential cutoff over the full analyzed energy 
range from 100 MeV to 820 GeV. For each of the three benchmark DGE models considered here, 
the power law index of the IGRB is ~ 2.3 and the cutoff energy is ~ 250 GeV (Table |4|. 

The total EGB is derived by adding resolved high-latitude LAT sources (taken to be primarily 
extragalactic) to the measured IGRB intensity. At an energy of 100 GeV, roughly half of the total 
EGB intensity has now been resolved into individual sources by the LAT, predominantly blazars of 
the BL Lacertae type. (The demographics of LAT sources detected at energies above 10 GeV are 


discussed in Ackermann et al. 2013c). The relative contribution of resolved sources becomes even 


more pronounced at energies exceeding 100 GeV. 

The intensities of the IGRB and the total EGB are compared to the first LAT measurement 
of the IGRB in Abdo et al. (2010b) in Figure [9j The two are compatible within the respective 
systematic uncertainties. Differences can be attributed to the combined effects of a more accurate 
estimate of the CR background at low energies, and changes in the Galactic foreground model. 
Importantly, the model for atmospheric secondaries has been refined to address discrepancies be- 
tween data and simulation. The revised background rate of misclassified CRs is up to 50% higher 
at a few hundred MeV than the older estimates. This change contributes to a reduced inte- 
grated IGRB intensity above 100 MeV of 7.2 ± 0.6 x 10“ 6 ph cnr 2 s -1 sr -1 in comparison to the 


1.03 ± 0.17 x 10 5 ph cm 2 s 1 sr 1 reported in Abdo et al. (2010b). 


The intensity resolved into individual sources at latitudes \b\ > 20° did not change substantially 


between the two measurements. This is consistent with the findings reported in Abdo et al. (2010c) 


and Ackermann et al. (2011b) that report a sky-averaged intensity of sources from Galactic latitudes 


\b\ > 10° of 4.0 x 10 -6 ph cm -2 s -1 sr -1 after one year and 4.4 x 10 -6 ph cm~ 2 s -1 sr -1 after 
two years of observations above 100 MeV. This difference corresponds to only ~ 5% of the IGRB 
intensity. 


Figure 10 places LAT measurements of the total EGB intensity in context with other mea- 
surements of the extragalactic X-ray and 7-ray backgrounds, together spanning nearly nine orders 
of magnitude in energy between 1 keV and 820 GeV. There is a good agreement between the to- 
tal EGB measured by the LAT and the previous measurement of the IGRB using EGRET data 
( Sreekumar et al.||1998 Strong et al.||2004 ) below 1 GeV. The IGRB measured by the LAT is lower 
than the EGRET IGRB measurement, as expected from the greatly superior sensitivity of the LAT 
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to resolve individual sources when compared to EGRET. 


As discussed in Section [lj numerous source populations as well as truly diffuse processes are 
expected to contribute to the EGB intensity. A detailed review of the expected contributions of 
specific source populations and diffuse processes is beyond the scope of this work. Instead, we focus 
on general constraints that can be applied to extragalactic 7 -ray source populations based on the 
EGB spectrum, taking into account the effects of EBL attenuation. Other efforts to statistically 


characterize the EGB properties considering the fluctuations of counts in spatial pixels (Malyshev 


& Hogg 2011) and two-point correlation functions (Ackermann et al. 2012a) have proven valuable 


for constraining the abundance of sources just below the LAT detection threshold, and similar 
techniques may be usefully applied to LAT data in the energy range > 100 GeV. 


In the interpretation that follows, we use the formalism outlined by Murase et al. (2007) 


and Inoue & Ioka (2012) to calculate both the EBL-attenuated primary signal as well as the 


electromagnetic cascade emission that would arise from a variety of generic cosmologically evolving 


source populations. We adopt the UV/optical/IR EBL model of Franceschini et al. (2008) based 
on observed galaxy counts, which is found to be consistent with spectral analyses of individually 
detected 7 -ray sources (e.g., Ackermann et al. 2012g). The populations are modeled as a collection 
of sources sharing the same intrinsic simple power-law spectral form with a common photon index, 
7 , and maximum energy, E max : 


dN/dE 


N 0 (E/E 0 ) ^ 1 E < E max 
0, E > E max 


(2) 


where E is the energy of photons emitted at the source. We model the evolving comoving volume 
emissivity (ph s _1 cm -3 ) of source populations without distinguishing between luminosity or density 
evolution of the sources. The emissivity evolution is either (1) parameterized by j oc (1 + z) 13 , or 
(2) follows the cosmic star formation rate density (Behroozi et al. 2013). We consider sources up to 
redshifts z = 2 and z = 3 in the two cases above, respectively. Our conclusions do not qualitatively 
change when using a lower maximum redshift oi z = 1. No attempt is made here to identify 
sources that could be individually resolved by the LAT versus those that blur into the unresolved 
background. 


The expected total EGB contributions corresponding to various source population scenarios 
are compared to the measured EGB spectrum > 20 GeV in Figures ED and [12] In each case, 
the contribution of the sources (i.e. , sum of primary and cascade components) to the total EGB 
has been normalized to match that of the measured total EGB intensity at 20 GeV. Figure 0 
shows source populations with parameterized comoving volume emissivities, while Figure \V2\ shows 
populations whose emissivity follows the comoving star formation rate density. In scenarios with 
photon index T = 1.5, the cascade component can dominate the primary component at energies 
< 100 GeV. 
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Fig. 9. — Comparison of the measured IGRB and total EGB intensities (foreground model A) to the 
first measurement of the IGRB in Abdo et al. (2010b) based on 10 months of LAT data. The error 
bars on the LAT measurements include the statistical uncertainty and systematic uncertainties 
from the effective area parametrization, as well as the CR background subtraction. Statistical and 
systematic uncertainties have been added in quadrature. The shaded bands indicate the systematic 
uncertainty arising from uncertainties in the Galactic foreground. The total EGB intensity is the 
sum of the IGRB and the intensity of the resolved LAT sources at high Galactic latitudes, |6| > 20°. 
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Fig. 10. — Comparison of the derived total EGB intensity (foreground model A) to other mea- 
surements of the X-ray and 7 -ray background. The error bars on the LAT measurement include 
the statistical uncertainty and systematic uncertainties from the effective area parametrization, as 
well as the CR background subtraction. Statistical and systematic uncertainties have been added 
in quadrature. The shaded band indicates the systematic uncertainty arising from uncertainties in 
the Galactic foreground. (Note that the EGRET measurements shown are measurements of the 
IGRB. However, EGRET was more than an order of magnitude less sensitive to resolve individual 
sources on the sky than the Fermi- LAT.) 
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Comoving Emissivity Density ~ (1+z) p , r = 1.5, E Max = 3 TeV Comoving Emissivity Density ~ (1+z) p , r = 1.5, E Mgx = 10 TeV 








Fig. 11. - EGB contributions of various source populations with comoving volume emissivity 
parameterized by j oc (1 + z)@ are compared to the measured EGB spectrum (foreground model 
A). Rows are differentiated by the assumed photon index of the sources. Left and right columns 
correspond to populations with maximum energies of 3 TeV and 10 TeV, respectively. Gray curves 
represent the intensity of primary 7 rays (attenuated by the EBL). Colored curves indicate the sum 
of the primary and cascade components. 


-40- 


Several patterns are apparent in Figures 11 and 12 First, source populations with negative 
evolution, especially those with very hard power-law spectra (r < 2) extending to multi- TeV 
energies, have difficulty fitting the high-energy break in the total EGB spectrum and are therefore 
unlikely to account for the EGB on their own. On the other hand, source populations whose 
emissivity evolves as the cosmic star-forming rate (corresponding to j3 ~ 3.25 at low redshifts, e.g., 
Hopkins &; Beacom 2006) also face challenges to match the shape of the total EGB spectrum alone, 
unless the spectral properties of those sources are finely tuned. The source populations which would 
most readily explain the measured total EGB spectrum from 100 MeV to nearly 1 TeV are those 
with photon indices matching that of the EGB below 100 GeV, namely ~ 2.3, and little or no 
evolution. In fact, the distribution of photon indices for individual LAT sources detected above 10 


GeV is also peaked near a value of 2.3 (Ackermann et al. 2013c). These and similar studies (e.g., 


Venters||2010 Inoue fc Ioka||20l2 Murase et al.|2012 ) demonstrate that the power-law shape of the 
total EGB spectrum with a single cutoff/break at ~ 250 GeV could in principle be explained by 
a single dominant extragalactic source population with relatively generic spectral properties and 
EBL attenuation. However, more sophisticated modeling efforts taking into account the specific 
properties of established extragalactic source classes are needed to fully understand how sources 
governed by such diverse physics produce a nearly featureless EGB spectrum over ~ 4 decades in 
energy. 

In addition to the implications for source classes comprising the EGB, a second aspect of this 
work is a further examination of the DGE (Section [4j). In the effort to accurately subtract the DGE 
and thereby isolate the fainter isotropic component, we considered a wider range of models for CR 
injection and propagation in the interstellar medium of the Milky Way than previously considered, 
e.g., in Ackermann et al. (2012b). None of the models tested here simultaneously satisfy constraints 


from both local CR measurements and observations of the high latitude 7 -ray sky, particularly in 
the case of IC emission (see Appendix |A|) . The lack of a clearly preferred DGE model is the largest 
single source of systematic uncertainty when measuring the IGRB intensity in the 100 MeV - 
100 GeV energy range with the LAT. Some of the modifications to commonly used CR injection 
and diffusion treatments investigated here may provide interesting future avenues of research. 


Further improvements over the Pass 7 event reconstruction and classification are required to 
extend EGB measurements with the LAT to both lower and higher energies. An extension to 
energies ~ 50 MeV may help to constrain the radiative processes contributing to the EGB, e.g., 
through the identification of the (redshifted) pionic spectral feature expected from the interactions 


at ~ 100 MeV. An extension to energies ~ 1 TeV would further clarify the spectra and evolution 
of sources that will be studied in detail with the extragalactic surveys of the High- Altitude Water 


of CR nuclei in all galaxies (e.g., Stecker & Venters 2011, Lacki et al. 2012; Chakraborty & Fields 


2013), although no indication of such a low-energy cutoff is present in the current measurement 
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Cherenkov observatory (HAWC^j and the Cherenkov Telescope Array (CTAj^j Both of these 
spectral extensions to the LAT IGRB measurement may be realized with future Pass 8 analyses 
(Atwood et al. [ [2013 ) . Additional insight regarding the sources of the EGB may come from ongoing 
studies of the extragalactic background of high-energy neutrinos (Aartsen et al. 2013), since the 
interactions of very-high and ultra-high energy CRs inevitably create fluxes of both 7 rays and 


neutrinos, the implications of which are discussed by, e.g., Ahlers et al. (2010), Berezinsky et al. 
(2011), Wang et al. (2011), Gelmini et al. (2012), and Murase et al. (2012). 


"http : / /www .hawc-observatory . org/ 
1( http : / /www . cta-observatory . org/ 
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Comoving Emissivity Density ~ SFR, E Max = 3TeV Comoving Emissivity Density ~ SFR, E Max = 10TeV 



Fig. 12. — EGB contributions of various source population scenarios with comoving volume emissiv- 
ity following the comoving star formation rate density are compared to the measured EGB spectrum 
(foreground model A). Each colored curve denotes a source population with a different assumed 
photon index. Left and right columns correspond to populations with maximum energies of 3 TeV 
and 10 TeV, respectively. Gray curves represent the intensity of primary 7 rays (attenuated by the 
EBL). Colored curves indicate the sum of the primary and cascade components. 
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A. Galactic foreground models 

This appendix summarizes the parameters used for modeling the diffuse Galactic foreground 
emission in the benchmark foreground models A, B and C. Please note that a customized version 
of GALPROP is needed to produce the models. The output files of the corresponding GALPROP 
runs are provided in an electronic data repository. A link to this repository is distributed in the 
electronic supplementary material. 


A.l. Foreground model A 


Foreground model A uses a parametrization of the distribution of pulsars in the Galaxy 


(Lorimer et al. 2006) as the distribution of CR sources (see Figure 13) where CR electrons and 
nuclei are injected into the interstellar medium. The diffusion coefficient for their propagation in 
the ISM is set to D xx = 7.0 X 10 ~ 28 m 2 s -1 (i?,/4 GV ) 0 ' 33 where R denotes the rigidity. The CRs 
are re-accelerated in the ISM with a re-acceleration strength (parametrized by the Alfven velocity 
va = 30 km s -1 ) that is constant throughout the Galaxy. The Galactic magnetic field is modeled 


according to Strong et al. (2011) using a local random field strength of 7.5 /rG. A CR halo size of 


5 kpc is chosen. All nuclei and electrons are injected with an energy spectrum that is a broken 
power law in rigidity. The power-law index of the injection spectrum of protons is 1.9 below 9 GV, 
2.45 between 9 GV and 240 GV, and 2.32 above 240 GV. For helium, we multiply the injection 
spectrum by R 01 , and the He/H fraction in the ISG was assumed to be 0.11. The second break 
in the energy spectrum and the harder injection spectrum for helium are motivated by the results 


of the measurements of the proton spectrum by the PAMELA satellite (Adriani et al. 2011b), 


and by the CR spectrum inferred from LAT observations of Earth limb 7 rays (Ackernrann et al 
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2014 


li 


The 7 -ray yield from the interactions of the CR with the ISG was calculated using the 


parametrization of Kamae et al. (2006) 


The power-law index of the injection spectrum of electrons is 1.5 below 5 GV, 2.85 between 
5 GV and 25 GV, and 2.32 above 25 GV. As in the case of nucleons, the breaks in the spectrum are 
introduced to obtain a good agreement of the model with measurements of the local CR electron 
spectrum ( Adriani et al.pOlla Abdo et al.||2009b JAckermann et al.||2012e ). The break at 5 GV is 


furthermore motivated by the shape of the Galactic synchrotron emission spectrum (Strong et al 


2011). Additionally, we assume a high-energy cutoff (implemented as a change in the injection 
index to 4 above 1.8 TeV) to remain in agreement with the H.E.S.S. measurements of the electron 


spectrum up to several TeV (Aharonian et al. 2008) 


The doubly broken injection spectra used in foreground model A improve the model/fit agree- 
ment in our maximum likelihood fit of the 7 -ray data. However, there is only negligible effect on the 


derived IGRB spectrum. Figure 14 shows a comparison between the expected 7 -ray spectra (from 
foreground model A) and the spectra obtained by the maximum likelihood procedure when fitting 
the model templates to the 7 -ray data (see also Section [3]). For the purpose of this comparison, 
we extend the upper bound of the energy range of the low-energy fit from 13 GeV to 51 GeV. 
The 7 rays arising from CR interactions with ISG and originating from the IC process are shown 
separately. The input model spectra are renormalized in the figure to allow for a better compar- 
ison of the predicted and the fitted spectral shapes. The renormalization factors are determined 
the energy band between 6.4 GeV and 51 GeV (cf. Section [3]). The numerical values of the 


m 


renormalization factors for the Hi + Hu and inverse Compton templates of our benchmark models 


are displayed in Figures 14, 16 and 18 


The predicted and observed spectra of 7 rays from CR interactions with ISG agree well at 
energies above a few GeV, besides a moderate renormalization factor. A harder spectrum is seen 
at energies below a few GeV in the fit compared to the model. In this energy range, the local 
interstellar spectrum of CRs is difficult to measure due to the effects of solar modulation. 


The IC emission model over-predicts the fitted IC emission at low energies by a factor of up 
to ~ 4, while at high energies it under-predicts the emission by a factor of 2.4. Again, the spectral 
shapes of the model and the fit to LAT data show good agreement at energies above a few GeV, so 
the (rescaled) model predictions from foreground model A can be used for the high-energy analysis, 
where the foreground model is fixed. 


11 We note that the AMS-02 collaboration has recently released preliminary data that does not confirm the hardening 
of the proton and helium spectra. A test showed that there are only negligible effects on the IGRB intensity derived 
in this work if we use a foreground model without a break at 240 GV in the nucleon spectrum. 
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A. 2. Foreground model B 


A significantly better agreement between the predicted and fitted IC emission is found for 
foreground model B, which includes an additional population of electron-only sources located near 
the Galactic center. In this model, the bulk of the IC emission arises from the electrons injected 
near the Galactic center, and the sources of CR nuclei do not produce a significant fraction of 
the local CR electron flux. The spatial distribution that is used for this additional electron-only 
source population is shown in Figure [l3j We further assume that the Galactic center sources inject 
electrons following a power-law spectrum with index 1 below 20 GV and 2.05 above 20 GV. We 
also assume the same high-energy cutoff as in foreground model A to maintain agreement with 
H.E.S.S. measurements (Aharonian et al. 2008). The additional electron source population must 
be located close to the Galactic center in order to produce bright enough IC emission to match 7 - 
ray observations without over-predicting the local CR electron spectrum at high energies (assuming 
that the model diffusion parameters are unchanged). 


Figure [T7T| shows the electron spectrum produced by model B in comparison to measurements, 
and the electron spectrum predicted by model A. There is a clear deficit of electrons below 20 GV in 
this model. However, these electrons could be easily supplied by a local source or source population 
of very soft electrons without a large impact on the total amount of IC emission. Natural candidates 
would be very old supernova remnants (SNRs) in the solar neighborhood, e.g., Loop I where we 
see high-energy 7 -ray emission from the shell in LAT data. It is, however, beyond the scope of this 
paper to speculate about and address the nature of such local soft electron sources. We ignore this 
potential local contribution for the IGRB analysis. 


The same propagation and injection parameters for nuclei are used in foreground model B as 
in model A. Modeled and fitted spectra are compared in Figure [16} The renormalization factor 
for 7 rays from the interactions of CRs with ISG increases from 1.5 in foreground model A to 1.7 
in foreground model B. This is a large change, and further investigations should be undertaken to 
understand if the model B value is still within the bounds of our current uncertainties concerning 
the total column density of the high-latitude gas, the 7 -ray emissivity of the ISG, and the gradient 
in the CR spectrum. For the purpose of this IGRB analysis, we accept the renormalized spectrum 
as a valid fit of the Galactic foreground. It can be further seen from Figure [TB] that the IC spectrum 
is now well described by the model both in normalization and in shape besides a small discrepancy 
at the lowest energies. 


A. 3. Foreground model C 

Foreground model C represents a class of models in which the CR diffusion and re-acceleration 
vary significantly throughout the Galaxy. Diffusion and re-acceleration are parametrized within 
the transport equation implemented in GALPROP via the spatial diffusion coefficient D xx and the 
Alfven velocity va- We use a simple model to vary the diffusion coefficient and Alfven velocity 
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by connecting their values to the strength of the regular and random Galactic magnetic fields. 
Following the approximation of Strong et al. (2007), the diffusion coefficient is set to 


D xx (R,r,z ) = D° xx 


( 8B{r , z) 

V sb q 


Be y 2 (R B q \ 
B(r,z)J \R 0 B(r,z) ) 


(Al) 


where B is the strength of the regular field and 8B corresponds to the strength of the random 
field at radius r from the Galactic center and height z above the Galactic plane. Bq and 5Bq 
denote their values at the position of the solar system. Rq = 4 GV is the reference rigidity, and 
D xx is the local diffusion coefficient at reference rigidity. The diffusion coefficient is constrained to 
D X x(R(h f, z) < 10 30 cm 2 s' 1 at reference rigidity. This constraint ensures that the mean free path 
of the CRs stays below the kpc scale for particles up to tens of TV in rigidit}p 2 | 


We parametrize the Alfven velocity as va oc B to t/y/~p, where B to t is the total magnetic field 
strength and p the density of ions in the interstellar medium. For the ion density, the same model 
is used as in the rest of this work (Gaensler et al. 2008). Simple models are assumed for the random 
and the regular magnetic field components with exponential scale heights and scale lengths. The 
regular magnetic field strength is assumed to be 4 pG at the position of the Sun, with a scale 
length of 11 kpc and a scale height of 4 kpc. The random field strength is assumed to be 4 pG, 
constant in the Galactic plane with a scale height of 2 kpc. The scale heights are in good agreement 
with scale heights found from equilibrium conditions (Kalberla & Kerp 1998). The field strengths 
are in qualitative agreement with recent studies of the Galactic synchrotron emission in Orlando] 
& Strong (2013), that take radio polarization into account. The extent of the CR halo for this 


model is set to 8 kpc; we note that constraints on the halo size found in earlier studies based on 
the 10 Be/ 9 Be ratio apply only to models with a static diffusion coefficient. Figure 17 shows the 
diffusion coefficient and the Alfven velocity as a function of the Galactocentric radius r and the 
height above the Galactic plane 2 . 


A customized version of the GALPROP code (see above) is used that allows the modeling of 
propagation scenarios in which the spatial and momentum diffusion are functions of radial distance 
from the Galactic center and height above the Galactic plane. As for foreground models A and 
B, GALPROP is used in its 2D mode that solves the transport equation on a 2D spatial grid in 
Galactocentric radius and height (r, z) around the Galactic center. The CR source distribution 
assumed in model C is more peaked toward the Galactic center than the pulsar distribution in 
model A (see Figure [13]). The high-energy injection spectra for CR electrons and protons are the 
same as for model A, power laws in rigidity with an index of 2.32. However, an injection spectrum 
with an index of 1.9 below 13 GV, and 2.40 between 13 GV and 240 GV is used for the CR protons 
in model C, slightly different from the spectrum used in model A. For the CR electrons the injection 


12 We tested that the specific choice of this upper bound on the diffusion coefficient is irrelevant by increasing 
the maximum value for D xx {Rq, r, z) by one order of magnitude. The effects on the predicted y-ray emission were 
negligible. 
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spectral index is 1.5 below 4.5 GV, and 2.70 between 4.5 GV and 25 GV. These modifications in 
the injection spectrum improve the agreement of the local CR spectra predicted by model C with 
measurements. 


Modeled and fitted 7 -ray spectra are compared in Figure 18 The renormalization factor for 
7 -rays from the interactions of CRs with the ISG is 1.5 as for model A. The total intensity of the 
IC emission is similarly under-predicted by a factor of 2.5 at energies above a few GeV, while it is 
over-predicted at low energies. However, an interesting aspect of this model is that it predicts a 
flatter CR gradient than models A and B, and therefore predicts a higher 7 -ray emissivity in the 
outer Galaxy. It was found in two other studies (Abdo et al. 2010 a, [Ackermann et al. 2011a) that 
the emissivity derived from LAT observations is indeed higher in the outer Galaxy than predicted 
by diffuse emission models of class A. We do not discuss the CR gradient of model C further here 
as it is not relevant for the IGRB analysis. 


-48- 



Fig. 13. — Parametrizations of the radial CR source distribution used in this work. The pulsar 


distribution is taken from Lorimer et al. (2006), the SNR distribution from Case & Bhattacharya 


(1998). The curves in the figure are normalized to unit integrated source density, the actual 


normalizations used in the models are derived from comparisons of predicted and measured local 
CR proton and electron intensities. 
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Energy [MeV] 

Fig. 14. — Fitted average intensity of the DGE in foreground model A for Galactic latitudes 
\b\ > 20°. Contributions from CR interactions with the ISG, and contributions from IC are shown 
separately. The normalizations of the two components are fitted individually in each energy bin. 
The GALPROP model spectrum that enters the fit is shown as dashed lines. The dashed lines are 
renormalized by the factors indicated in the legend. 
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Fig. 15. — Predicted local CR electron spectrum of foreground models A and B in comparison to 
measurements of the spectrum. LAT data are taken from Ackermann et al. (2010), PAMELA data 


from Adriani et al. (2011a), AMS-01 data from Aguilar et al. (2002), H.E.S.S. data from Aharonian 


et al. (2009) and Egberts et al. (2011). 
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Fig. 16. 
tion. 


Fitted average intensity of the DGE in foreground model B. See Figure 14 for a descrip- 



r [kpc] r [kpc] 


Fig. 17. — Diffusion coefficient (left) and Alfven velocity (right) used in foreground model C. Both 
are functions of distance from the Galactic center (r) and height above the Galactic plane (z). 


Alfven speed v [ km s' 



E 2 dN/dE [cm 2 s 1 sr 1 MeV] 
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Fig. 18. — Fitted average intensity of the DGE in foreground model C. See Figure [14] for a descrip- 
tion. 
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B. Residual maps 


Figure [19] shows residual maps of the relative deviations in the number of expected and observed 
counts when using foreground model A for the DGE. The first map shows the residual for all counts 
above 100 MeV, the second map shows the residual for counts above 13 GeV. Multiple structures 
are visible in the former map while the latter is dominated by the Fermi Bubbles. 


Figure 20 visualizes the difference in the predicted 7 -ray emission when using Galactic fore- 
ground models B or C instead of model A. These differences are more prominent at energies above 
a few GeV, where the IC emission contributes a larger fraction of the total 7 -ray emission than at 
few hundreds of MeV. Therefore, the relative deviation in predicted counts above 3 GeV is shown 
in the maps when using foreground models B and C in the fit with respect to using foreground 
model A. 


For foreground model B, a higher 7 -ray intensity is predicted close to the Galactic center region 
arising from the IC emission of electrons that originate from the additional source population we 
introduced in this model. For foreground model C, regions with higher intensity can be observed 
towards the outer Galaxy, reflecting the more efficient transport of CRs into the outer Galaxy by 
the modified propagation scheme used in model C. 


-54- 



180.0 


120.0 


60.0 0.0 300.0 

Galactic longitude [deg] 


240.0 


180.0 


( Counts - Model ) / Model 


- 0.2 - 0.15 - 0.1 - 0.05 0 0.05 0.1 0.15 0.2 



180.0 


120.0 


60.0 0.0 300.0 

Galactic longitude [deg] 


240.0 


180.0 


( Counts - Model ) / Model 


- 0.8 - 0.6 - 0.4 - 0.2 0 0.2 0.4 0.6 0.8 


Fig. 19. — Residual maps for foreground model A used in this analysis. The fractional difference 
in counts between the actual data and the fitted model is shown in the figures. Upper: All counts 
above 100 MeV are included in the map. The pixel size is 0.8 deg 2 (HEALPix order 6). Lower: 
Only counts above 13 GeV are included in the map. The pixel size has been increased to 13 deg 2 
(HEALPix order 4) to account for the reduced count statistics at higher energies. 
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Fig. 20.- Differences in foreground models A, B, and C. The fractional difference in the number 
of predicted counts above 3 GeV between the models using the alternative Galactic foregrounds, 
B and C and the model using Galactic foreground A are displayed in the maps. The pixel size is 
0.8 deg 2 (HEALPix order 6). Upper: Fractional difference for model B. Lower: Fractional difference 
for model C. 
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C. Criteria for masking regions in the likelihood fit 


Approximately 90% of the ISG is atomic (Hi), ionized (Hu), or molecular hydrogen gas (H 2 ). 
The remainder is mostly helium, and is generally assumed to be mixed uniformly with the hydrogen. 
The column density and distribution of atomic gas in the Galaxy can be estimated from surveys 


of the 21 cm hyperfine-structure transition line of the hydrogen atom (Kalberla et al. 2005). The 
distribution of molecular hydrogen gas can be inferred indirectly from surveys of the 2.6 mm J(l— >-0) 
transition of the CO molecule by assuming a proportionality (usually called the Aco factor) between 


the intensity of this line integrated over frequency (Rco) and the column density of H 2 gas (Dame 
et al. 2001). In both cases, the velocity component of the gas parallel to the line of sight is measured 
via the doppler shift of the transition line. When combined with a model of the Galaxy rotation 
curve, the observed velocity can be converted into a measurement of Galactocentric radius, allowing 
for a determination of the gas density distribution along the line of sight. In this work, we use 


the gas distribution model of (Ackermann et al. 2012b), where the total gas column density is 
distributed in 17 Galactocentric rings spanning the radial range from 0 to 50 kpc. Due to the 
small scale heights of the gas (a few tens of pc for H 2 and a few hundreds of pc for Hi), most of 
the gas outside of our local Galactic neighborhood will appear concentrated around the Galactic 
plane. Even in our local neighborhood, most of the H 2 gas is concentrated in isolated clouds at low 
Galactic latitudes. 

We use this fact to exclude regions of the sky from our likelihood fit that have a significant 
column density of H 2 gas (Wqo > 2.5 K km s -1 ) along the respective line of sight as well as lines 
of sight with a significant column density of Hi gas (IVhi > 5 x 10 20 cnr 2 ) located beyond our local 


solar neighborhood (8 kpc < r <10 kpc) according to the gas distribution model in Ackermann 


et al. (2012b). Independent of the gas column densities found, Galactic latitudes |6| <5° are also 
excluded. The exclusion of such regions simplifies the likelihood analysis considerably. Specifically, 
models of the DGE in the remaining mostly high-latitude parts of the sky do not depend on 
assumptions about the CO-to-H 2 conversion factor Xqo, or on how accurately we model variations 
of the CR density throughout the Galaxy leading to variations of the Hi emissivity for gas at larger 
distances. The excluded regions cover a total of 17% of the sky. 


REFERENCES 

Aartsen, M. G., Abbasi, R., Abdou, Y., et al. 2013, Science, 342, 1242856 
Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009a, Phys. Rev. D, 80, 122004 
— . 2009b, Physical Review Letters, 102, 181101 
— . 2009c, Astroparticle Physics, 32, 193 
— . 2010a, ApJ, 710, 133 


- 57 - 


— . 2010b, Physical Review Letters, 104, 101101 
— . 2010c, ApJ, 720, 435 
— . 2011, ApJ, 734, 116 
— . 2012, ApJ, 758, 140 

Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17 

Ackermann, M., Johannesson, G., Digel, S. W., et al. 2009, AIP Conf. Proc., 1085, 763 

Ackermann, M., Ajello, M., Atwood, W. B., et al. 2010, Phys. Rev. D, 82, 092004 

Ackermann, M., Ajello, M., Baldini, L., et al. 2011a, ApJ, 726, 81 

Ackermann, M., Ajello, M., Allafort, A., et al. 2011b, ApJ, 743, 171 

Ackermann, M., Ajello, M., Albert, A., et al. 2012a, PhRvD, 85, 083007 

Ackermann, M., Ajello, M., Atwood, W. B., et al. 2012b, ApJ, 750, 3 

Ackermann, M., Ajello, M., Allafort, A., et al. 2012c, ApJ, 755, 164 

— . 2012d, Astroparticle Physics, 35, 346 

— . 2012e, Physical Review Letters, 108, 011103 

Ackermann, M., Ajello, M., Albert, A., et al. 2012f, ApJS, 203, 4 

Ackermann, M., Ajello, M., Allafort, A., et al. 2012g, Science, 338, 1190 

— . 2013a, ApJ, 765, 54 

— . 2013b, ApJ, 778, 82 

— . 2013c, ApJS, 209, 34 

Ackermann, M., Ajello, M., Albert, A., et al. 2014, ArXiv e-prints, astro-ph: 1403.5372 

Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2011a, Physical Review Letters, 106, 
201101 

— . 2011b, Science, 332, 69 

Aguilar, M., Alcaraz, J., Allaby, J., et al. 2002, Phys. Rep., 366, 331 

Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2008, Physical Review Letters, 
101, 261104 



- 58 - 


Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009, A&A, 508, 561 

Ahlers, M., Anchordoqui, L. A., Gonzalez-Garcia, M. C., Halzen, F., & Sarkar, S. 2010, APh, 34, 
106 

Ajello, M., Greiner, J., Sato, G., et al. 2008, ApJ, 689, 666 

Ajello, M., Shaw, M. S., Romani, R. W., et al. 2012, ApJ, 751, 108 

Ajello, M., Romani, R. W., Gasparrini, D., et al. 2014, ApJ, 780, 73 

Atwood, W., Albert, A., Baldini, L., et al. 2013, ArXiv e-prints, astro-ph: 1303.3514 

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071 

Bechtol, K. 2012, PhD thesis, Stanford University 

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 

Berezinskii, V. S., & Smirnov, A. I. 1975, Ap&SS, 32, 461 

Berezinsky, V., Gazizov, A., Kachelriefi, M., & Ostapchenko, S. 2011, PhLB, 695, 13 

Bergstrom, L., Edsjo, J., &; Ullio, P. 2001, PhRvL, 87, 251301 

Bregeon, J., Charles, E., & Wood, M. 2013, ArXiv e-prints, astro-ph: 1304.5456 

Casandjian, J.-M., & Grenier, I. 2009, ArXiv e-prints, astro-ph:0912.3478 

Case, G. L., & Bhattacharya, D. 1998, ApJ, 504, 761 

Chakraborty, N., & Fields, B. D. 2013, ApJ, 773, 104 

Churazov, E., Sunyaev, R., Revnivtsev, M., et al. 2007, A&A, 467, 529 

Clark, G. W., Garmire, G. P., &: Kraushaar, W. L. 1968, ApJ, 153, L203 

Colafrancesco, S., &: Blasi, P. 1998, APh, 9, 227 

Coppi, P. S., & Aharonian, F. A. 1997, ApJ, 487, L9 

Dame, T. M., Hartmann, D., Sz Thaddeus, P. 2001, ApJ, 547, 792 

de Palma, F., Brandt, T. J., Johannesson, G., & Tibaldo, L. 2013, ArXiv e-prints, astro- 
ph: 1304. 1395 

Dernier, C. D. 2007, AIP Conf. Proc., 921, 122 

Di Mauro, M., Calore, F., Donato, F., Ajello, M., & Latronico, L. 2014, ApJ, 780, 161 



- 59 - 


Dobler, G., Finkbeiner, D. P., Cholis, I., Slatyer, T., &; Weiner, N. 2010, ApJ, 717, 825 
Egberts, K., et al. 2011, Nuclear Instruments and Methods in Physics Research A, 630, 36 
Faucher-Giguere, C.-A., & Loeb, A. 2010, JCAP, 1, 5 
Feldmann, R., Hooper, D., & Gnedin, N. Y. 2013, ApJ, 763, 21 
Fermi-LAT Collaboration. 2014, submitted to ApJ 

Fichtel, C. E., Hartman, R. C., Kniffen, D. A., et al. 1975, ApJ, 198, 163 

Fichtel, C. E., Simpson, G. A., & Thompson, D. J. 1978, ApJ, 222, 833 

Fields, B. D., Pavlidou, V., & Prodanovic, T. 2010, ApJ, 722, L199 

Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837 

Fukada, Y., Hayakawa, S., Kasahara, I., et al. 1975, Nature, 254, 398 

Gaensler, B. M., Madsen, G. J., Chatterjee, S., & Mao, S. A. 2008, PASA, 25, 184 

Gelmini, G. B., Kalashev, O., & Semikoz, D. V. 2012, J. Cosmology Astropart. Phys., 1, 44 

Gendreau, K. C., Mushotzky, R., Fabian, A. C., et al. 1995, PASJ, 47, L5 

Gorski, et al. 2005, ApJ, 622, 759 

Grenier, I. A., Casandjian, J.-M., & Terrier, R. 2005, Science, 307, 1292 
Gruber, D. E., Matteson, J. L., Peterson, L. E., & Jung, G. V. 1999, ApJ, 520, 124 
Haslarn, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E. 1982, A&AS, 47, 1 
Healey, S. E., Romani, R. W., Taylor, G. B., et al. 2007, ApJS, 171, 61 
Hopkins, A. M., & Beacom, J. F. 2006, ApJ, 651, 142 
Inoue, Y. 2011, ApJ, 733, 66 

Inoue, Y., & Ioka, K. 2012, Phys. Rev. D, 86, 023003 
James, F., & Roos, M. 1975, Cornput . P hys . Coirimun 10, 343 
Johannesson, G., &: Orlando, E. 2013, ArXiv e-prints, astro-ph: 1307.0197 
Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775 
Kalberla, P. M. W., & Kerp, J. 1998, A&A, 339, 745 

Karnae, T., Karlsson, N., Mizuno, T., Abe, T., & Koi, T. 2006, Astrophys. J., 647, 692 



- 60 - 


Keshet, U., Waxman, E., & Loeb, A. 2004, Journal of Cosmology and Astro-Particle Physics, 4, 6 
Kinzer, R. L., Jung, G. V., Gruber, D. E., Matteson, J. L., & Peterson, L. E. 1997, ApJ, 475, 361 
Kraushaar, W. L., Clark, G. W., Garmire, G. P., et al. 1972, ApJ, 177, 341 

Kulkarni, S. R., & Heiles, C. 1988, Neutral hydrogen and the diffuse interstellar medium, ed. K. I. 
Kellermann &; G. L. Verschuur, 95-153 

Lacki, B. C., Horiuchi, S., Sz Beacom, J. F. 2012, ArXiv e-prints, astro-ph: 1206.0772 
Loeb, A., & Waxman, E. 2000, Nature, 405, 156 

Lorimer, D. R., Faulkner, A. J., Lyne, A. G., et al. 2006, MNRAS, 372, 777 
Makiya, R., Totani, T., & Kobayashi, M. A. R. 2011, ApJ, 728, 158 
Malyshev, D., Sz Hogg, D. W. 2011, ApJ, 738, 181 
Moskalenko, I. V., Sz Porter, T. A. 2009, ApJ, 692, L54 
Moskalenko, I. V., Porter, T. A., Sz Digel, S. W. 2006, ApJ, 652, L65 
Murase, K., Asano, K., Sz Nagataki, S. 2007, ApJ, 671, 1886 

Murase, K., Beacom, J. F., Sz Takami, H. 2012, J. Cosmology Astropart. Phys., 8, 30 
Nolan, P. L., Abdo, A. A., Ackermann, M., et al. 2012, ApJS, 199, 31 
Orlando, E., Sz Strong, A. 2013, MNRAS, 436, 2127 
Orlando, E., Sz Strong, A. W. 2007, Ap&SS, 309, 359 
— . 2008, A&A, 480, 847 

Pavlidou, V., Sz Fields, B. D. 2002, ApJ, 575, L5 

Porter, T. A., Moskalenko, I. V., Strong, A. W., Orlando, E., Sz Bouchet, L. 2008, ApJ, 682, 400 

Revnivtsev, M., Gilfanov, M., Sunyaev, R., Jahoda, K., Sz Markwardt, C. 2003, ASzA, 411, 329 

Schlegel, D. J., Finkbeiner, D. P., Sz Davis, M. 1998, ApJ, 500, 525 

Singal, J., Petrosian, V., Sz Ajello, M. 2012, ApJ, 753, 45 

Sreekumar, P., Bertsch, D. L., Dingus, B. L., et al. 1998, ApJ, 494, 523 

Stecker, F. W., Sz Venters, T. M. 2011, ApJ, 736, 40 

Strong, A. W., Moskalenko, I. V., Sz Ptuskin, V. S. 2007, Annual Review of Nuclear and Particle 
Science, 57, 285 



61 - 


Strong, A. W., Moskalenko, I. V., & Reimer, 0. 2000, Astrophys. J., 537, 763 
— . 2004, Astrophys. J., 613, 956 

Strong, A. W., Orlando, E., & Jaffe, T. R. 2011, A&A, 534, A54 

Su, M., Slatyer, T. R., & Finkbeiner, D. P. 2010, ApJ, 724, 1044 

Taylor, J. E., & Silk, J. 2003, MNRAS, 339, 505 

Thompson, T. A., Quataert, E., & Waxrnan, E. 2007, ApJ, 654, 219 

Ullio, P., Bergstrom, L., Edsjo, J., & Lacey, C. 2002, PhRvD, 66, 123502 

Venters, T. M. 2010, ApJ, 710, 1530 

Vladimirov, A. E., Digel, S. W., Johannesson, G., et al. 2011, Computer Physics Communications, 
182, 1156 

Wang, X.-Y., Liu, R.-Y., & Aharonian, F. 2011, ApJ, 736, 112 

Watanabe, K., Hartmann, D. H., Leising, M. D., et al. 1997, in American Institute of Physics 
Conference Series, Vol. 410, Proceedings of the Fourth Compton Symposium, ed. C. D. 
Dernier, M. S. Strickman, & J. D. Kurfess, 1223-1227 

Weidenspointner, G., Varendorff, M., Kappadath, S. C., et al. 2000, in American Institute of 
Physics Conference Series, Vol. 510, American Institute of Physics Conference Series, ed. 
M. L. McConnell & J. M. Ryan, 467-470 

Wolleben, M. 2007, ApJ, 664, 349 


This preprint was prepared with the AAS LATgX macros v5.2. 



